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ABSTRACT 


The  numerical  solution  of  electromagnetic  scattering  problems  involves 
approximating  an  exact  equation  by  a  finite-dimensional  matrix  equation.  The 
use  of  an  iterative  algorithm  to  solve  the  matrix  equation  sometimes  results  in 
a  considerable  savings  in  computer  memory  requirements.  For  a  fixed  amount  of 
computer  memory,  this  approach  permits  the  analysis  of  scatterers  that  are  an 
order  of  magnitude  larger  electrically. 

Iterative  algorithms  of  the  conjugate  gradient  class  are  examined  and 
applied  to  a  variety  of  typical  electromagnetic  scattering  problems,  in  order  to 
evaluate  their  performance  in  practice.  In  contrast  with  the  simple  iterative 
algorithms  used  in  the  past,  which  often  diverged  when  applied  to  electromag¬ 
netics  problems,  these  algorithms  never  diverge  and  usually  converge  at  a  quick 
rate . 

Depending  on  the  geometry  of  the  scatterer  under  consideration,  it  may  be 
possible  to  build  symmetries  into  the  matrix  representation  and  effect  the 
necessary  storage  reduction.  Two  distinct  approaches  for  creating  these  sym¬ 
metries  are  examined.  An  alternate  procedure,  which  requires  some  of  the  matrix 
elements  to  be  regenerated  as  needed  by  the  iterative  algorithm  in  use,  does  not 
rely  on  symmetries  and  is  applicable  to  a  larger  set  of  geometries.  Both  proce¬ 
dures  are  applied  to  several  scattering  problems.  Execution  time  comparisons 
show  that  the  approaches  based  on  symmetries  are  the  most  efficient,  and  that 
both  procedures  can  be  superior  to  noniterative  techniques  for  large  scatterers. 
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1.  INTRODUCTION 

Electromagnetic  scattering  problems  comprise  a  variety  of  physical 
phenomena,  from  the  interaction  of  radio  waves  with  aircraft  to  the  biological 
effects  of  microwaves  on  the  human  body.  Three  distinct  approaches  have  evolved 
for  the  analysis  of  scattering  problems.  The  first,  exact  analytical  techniques, 
include  approaches  such  as  separation  of  variables  [1],  [2]  and  the  Wiener-Hopf 
method  [3],  Although  exact  techniques  have  been  developed  and  applied  to 
electromagnetics  since  the  appearance  of  Maxwell's  work  [4],  they  remain  limited 
to  problems  which  can  be  described  by  relatively  simple  geometries.  The  second, 
approximate  analytical  methods,  such  as  the  variational  methods  [5]  and  the 
asymptotic  approaches  including  geometric  optics  [6]  and  geometric  theory  of 
diffraction  [7],  can  sucessfully  treat  many  types  of  problems  not  amenable  to 
exact  solution.  The  asymptotic  methods  are  usually  limited  to  the  analysis  of 
electrically  large  scatterers  whose  geometry  can  be  described  in  terms  of  the 
few  canonical  shapes  for  which  diffraction  coefficients  are  available.  The  suc¬ 
cess  of  these  methods  is  usually  highly  dependent  on  the  intuition  and  skill  of 
the  user,  and  systematic  estimates  of  the  accuracy  are  usually  impossible  to 
obtain.  The  final  category,  numerical  solutions,  encompasses  the  many  computer- 
aided  approaches  developed  over  the  past  several  decades  [8],  [9].  Numerical 
solutions  are  not  fundamentally  restricted  to  scatterers  with  certain  canonical 
shapes  or  materials,  and  in  principle  they  can  be  carried  out  to  obtain  any 
level  of  accuracy.  They  are,  however,  limited  in  practice  by  the  available  com¬ 
puter  resources.  Because  of  this  limitation,  conventional  numerical  methods  are 
best  suited  for  the  analysis  of  electrically  small  scatterers,  i.e.,  scatterers 
whose  maximum  dimensions  are  several  wavelengths  or  less.  Recently,  specialized 
numerical  techniques  have  been  introduced  for  the  efficient  treatment  of  larger 
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scatterers  [10]  -  [13].  The  purpose  of  this  thesis  is  to  review  these  recent 
developments  and  extend  these  techniques  to  other  applications. 

Every  numerical  solution  process  can  be  separated  into  three  parts.  First, 
the  physical  problem  must  be  cast  into  the  form  of  a  mathematical  equation. 

Since  few  equations  arising  in  practice  are  convenient  to  solve  analytically, 
the  second  part  of  the  process  involves  replacing  the  original  equation  by  a 
finite  dimensional  matrix  equation,  which  is  then  amenable  to  solution.  This  is 
known  as  discretization.  The  final  part  of  the  process  is  to  solve  the  matrix 
equation,  either  by  direct  methods  such  as  Gaussian  elimination  or  by  iterative 
methods . 

Initially,  the  problem  may  be  posed  as  a  differential  equation  with  boun¬ 
dary  conditions,  an  integral  equation,  or  some  variational  principle  to  be 
satisfied.  Furthermore,  the  problem  may  be  posed  in  the  time  or  frequency 
domain.  The  most  popular  approach  for  the  numerical  treatment  of  electromagnetic 
scattering  problems  has  been  the  integral  equation  formulation  [14]  -  [16], 
although  progress  has  been  reported  using  other  formulations  [17],  Often,  these 
integral  equations  are  convolutional  in  form,  and  this  property  is  of  central 
importance  for  many  of  the  specialized  methods  to  be  considered  below. 

The  discretization  process  may  involve  the  direct  sampling  of  pertinent 
quantities  over  a  grid  of  points,  often  used  in  connection  with  finite- 
difference  approximations  to  derivatives  [18].  If  the  problem  is  posed  in  terms 
of  a  variational  principle,  the  discretization  may  take  the  form  of  the  Ritz 
procedure  where  the  unknown  quantity  is  approximated  by  several  trial  functions 
[19],  [20].  The  approach  widely  used  to  discretize  the  integral  equations  of 
electromagnetics  is  a  generalization  of  the  Ritz  procedure  known  as  the 
method  of  moments  (MoM)  [21]  or  the  weighted  residual  method  [22].  The  MoM 
received  widespread  attention  during  the  1960's,  and  excellent  reviews  of  its 
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early  use  in  electromagnetics  are  provided  by  Richmond  [23],  Harrington 
[24]  -  [26],  and  Tanner  and  Andreasen  [27].  The  MoM  procedure  requires  the 
unknown  function  to  be  approximated  by  N  expansion  functions,  and  reduces  the 
task  to  that  of  finding  N  unknown  coefficients.  The  resulting  equation  is  made 
orthogonal  to  N  testing  functions,  to  yield  an  NxN  matrix  equation.  This  proce¬ 
dure  is  illustrated  in  Chapter  5,  and  additional  information  regarding  the  MoM 
may  be  found  in  several  recent  textbooks  [8],  [25],  [28],  [29]. 

Algorithms  for  the  solution  of  matrix  equations  are  also  available  in  many 
texts  [30]  -  [34].  Conventional  numerical  analysis  usually  incorporates  a 
direct  method,  such  as  Gaussian  elimination,  to  solve  the  matrix  equation.  The 
specialized  techniques  developed  to  treat  large  scatterers  usually  use  an  itera¬ 
tive  method  to  solve  the  system  because  general  purpose  elimination  algorithms 
require  the  NxN  matrix  to  be  stored  in  computer  memory.  Computer  memory 
consists  of  "fast  access"  in-core  memory  and  "slow  access"  out-of-core 
memory,  the  latter  including  peripheral  devices  such  as  disk  and  magnetic 
tape  facilities.  Because  of  physical  limitations,  only  relatively  small  order 
(i.e.,  orders  of  several  hundred)  matrices  can  be  stored  in  most  computers 
without  some  use  of  out-of-core  memory.  Unfortunately,  the  order  of  the  matrix 
increases  with  the  electrical  size  of  the  scatterer,  and  frequently  exceeds  the 
fast-access  limits  of  a  given  machine.  Although  advances  in  computer  architec¬ 
ture  continuously  improve  the  efficiency  of  transferring  data  to  and  from 
peripheral  storage  facilities,  such  transfers  remain  extremely  slow  compared 
to  in-core  access  times.  Therefore,  the  specialized  techniques  of  interest 
attempt  to  reduce  the  necessary  storage  to  a  level  which  can  be  handled  by  in- 
core  computer  memory  alone,  or  at  least  minimize  the  transfers  to  out-of-core 
devices.  Iterative  algorithms  do  not  require  the  NxN  matrix  to  be  stored  in  com¬ 
puter  memory,  and  thus,  are  highly  compatible  with  the  specialized  techniques 
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because  they  can  take  into  account  any  sparseness  or  redundancy  of  the  matrix  to 
reduce  the  storage  burden. 

Originally,  specialized  techniques  for  large  scatterers  incorporated 
"simple"  iterative  algorithms  of  the  Jacobi,  Gauss-Seidel ,  or  Successive 
Over-Relaxation  (SOR)  variety,  which  are  described  in  detail  in  the  literature 
[30]  -  [34].  These  algorithms  can  be  problematic,  because  they  sometimes 
diverge.  More  "sophisticated"  iterative  algorithms,  such  as  the  gradient 
methods  [30],  [33],  are  more  complicated  to  implement  but  never  diverge.  This 
report  will  only  consider  three  algorithms  based  on  gradient  methods,  all  of 
which  are  developed  in  Chapter  2. 

There  is  one  major  drawback  to  the  use  of  iteration  for  computational 
electromagnetics.  The  efficiency  of  iteration  is  poor  compared  to  direct  methods 
of. solution  whenever  multiple  systems  involving  the  same  NxN  matrix  must  be 
treated.  Because  direct  methods  in  effect  generate  an  inverse  matrix,  they  can 
efficiently  solve  a  matrix  equation  involving  any  number  of  right-hand  sides.  To 
date,  no  general,  systematic  iterative  procedure  is  available  with  similar  effi¬ 
ciency.  The  relative  inefficiency  of  iteration  holds  to  some  degree  in  spite  of 
the  fact  that  out-of-core  storage  procedures  may  be  necessary  for  the  implemen¬ 
tation  of  a  direct  method.  For  electromagnetic  scattering  problems,  there  is  a 
trade-off  between  these  two  approaches  and  this  trade-off  apparently  has  not 
been  studied  in  detail  to  date. 

All  of  the  well-known  frequency-domain  based  techniques  for  the  efficient 
analysis  of  electrically  large  scatterers  have  incorporated  three  features 
[101  ■  [13].  First,  they  involve  a  convolutional  integral  equation  formulation. 
Second,  they  are  restricted  to  geometries  which  can  be  discretized  in  such  a 
manner  as  to  preserve  the  convolution  form  in  the  matrix  equation.  (Such  a 
system  consists  of  one  or  more  discrete  convolution  operations,  meaning  that  it 
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possesses  a  significant  degree  of  redundancy.)  Finally,  they  incorporate  itera¬ 
tive  algorithms  to  solve  the  matrix,  and  attempt  to  implement  these  algorithms 
so  as  to  minimize  the  computer  memory  requirements. 

Although  this  report  will  not  discuss  time-domain  based  techniques  in  any 
detail,  they  deserve  some  mention  because  they  have  also  been  applied  to  the 
efficient  analysis  of  electrically  large  scatterers  [35)  -  [38].  Time-domain 
approaches  can  require  considerably  less  computer  memory  and  execution  time  than 
conventional  frequency-domain  based  methods,  and  appear  to  be  comparable  in 
efficiency  to  the  specialized  frequency-domain  methods  considered  here.  General 
reviews  of  numerical  time-domain  techniques  are  provided  by  Mittra  [39],  Bennett 
and  Ross  [40],  and  Miller  and  Landt  [41]. 

The  first  specialized  frequency-domain  based  technique  for  the  analysis  of 
electrically  large  scatterers  was  published  in  a  comprehensive  report  by 
Bojarski  in  1971  [10].  Although  Bojarski's  "K-space"  formulation  made  extensive 
mention  of  the  Fourier  transform  and  the  fast-Fourier  transform  (FFT)  algorithm, 
the  distinguishing  features  of  the  approach  as  a  numerical  technique  included 
its  use  of  a  convolutional  integral  equation,  an  evenly-spaced  sampling  grid 
which  preserved  the  convolution  form  in  the  matrix  equation,  and  a  "simple" 
iterative  algorithm  to  solve  the  system.  Bojarski  made  extensive  use  of  the  FFT, 
which  allows  a  large  order  discrete  convolution  to  be  performed  more  efficiently 
than  possible  by  direct  matrix  multiplication  [42],  [43].  However,  the  FFT  can 
only  be  directly  applied  to  an  infinite-periodic  seouence  of  numbers;  the  treat¬ 
ment  of  non-periodic  sequences  requires  the  incorporation  of  zero-padding  [42]. 
Because  of  this,  any  use  of  the  FFT  algorithm  for  non-periodic  scattering 
problems  increases  the  array  sizes  and  the  computer  memory  requirements 
cons  iderab ly . 
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Ko  and  Mittra  [44],  and  later  Kastner  and  Mittra  [11],  [45],  [46],  developed 
a  technique  similar  to  Bojarski's  that  they  applied  to  a  variety  of  applica¬ 
tions.  Entitled  the  "Spectral-Iterative"  technique  (SIT),  it  differed  from  the 
"K-space"  method  in  two  respects.  First,  the  "K-space"  method  required  a  three- 
dimensional  FFT  and  three-dimensional  zero-padding  to  treat  a  three-dimensional 
non-periodic  scatterer.  The  SIT  required  only  a  two-dimensional  FFT  to  treat  a 
three-dimensional  geometry,  and  thus  eliminated  some  of  the  cumbersome  zero¬ 
padding.  Since  the  SIT  required  the  discrete  convolution  to  be  summed  explicitly 
in  the  remaining  dimension,  it  would  require  somewhat  more  computation  per 
iteration  step  than  the  "K-space"  method  for  large  scatterers.  It  is  important 
to  realize  that  an  explicit  summation  can  be  substituted  for  the  FFT  whenever 
desired,  and  may  be  necessary  if  the  additional  storage  constraints  imposed  by 
the  FFT  algorithm  exceed  the  fast-access  limits  of  a  given  machine.  The  second 
difference  between  the  SIT  and  the  "K-space"  method  lies  in  the  discretization 
employed.  The  discretization  used  within  the  SIT  required  the  direct  sampling  of 
the  analytical  Fourier  transform  of  the  kernel  (Green's  function)  appearing  in 
the  integral  equation,  as  opposed  to  a  sampling  of  the  kernel  itself  as  was  done 
in  the  "K-space"  method.  The  difference  between  these  two  approaches  is  mani¬ 
fested  in  the  difference  between  the  Fourier  transform  and  the  FFT,  although 
other  factors  are  also  involved. 

Since  the  FFT  treats  a  sequence  as  periodic  in  both  the  original  and  the 
transform  domains,  the  validity  of  mixing  the  analytical  transform  and  the  FFT 
for  non-periodic  functions  comes  under  question.  Because  of  this,  Bojarski  spe¬ 
cifically  recommended  against  the  type  of  discretization  that  was  later  used 
with  the  SIT  [10].  In  a  variety  of  applications,  some  of  which  are  illustrated 
in  Chapter  5,  it  is  desirable  to  sample  the  Fourier  transform  directly  because 
the  Green's  function  is  not  well  suited  for  numerical  computation.  In  order  to 
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ascertain  the  validity  of  the  type  of  discretization  in  which  this  is  done,  and 
to  develop  guidelines  for  its  use.  Chapter  5  presents  a  detailed  comparison  of 
the  two  discretization  procedures. 

Although  both  the  "K-space"  method  and  the  SIT  were  applied  to  a  variety  of 
problems,  both  methods  were  inconvenient  to  use  in  practice  because  the 
"simple"  iterative  algorithm  used  to  solve  the  matrix  equation  often  diverged. 
Recently,  Sultan  and  Mittra  incorporated  a  "sophisticated"  iterative  algorithm, 
the  conjugate  gradient  method  (CGM) ,  into  the  SIT  procedure  [13].  The  CGM  has 
only  recently  received  widespread  attention  in  electromagnetics,  and  thus  its 
performance  in  practice  is  largely  undocumented.  In  an  attempt  to  alleviate  this 
situation.  Chapter  2  describes  the  CGM  and  two  related  algorithms  in  detail,  and 
Chapter  3  discusses  the  performance  of  the  algorithms  for  a  variety  of  scat¬ 
tering  problems. 

Other  iterative-based  methods  have  been  developed  to  efficiently  treat 
scattering  problems.  Nyo  and  Harrington  [12],  [47]  and  Borup  and  Gandhi  [48] 
applied  the  MoM  discretization  to  geometries  which  preserved  the  discrete  con¬ 
volutional  form  of  the  original  equation.  The  resulting  matrix  equations  were 
solved  with  a  simple  iterative  algorithm,  although  Borup  and  Gandhi  later  incor¬ 
porated  the  CGM  into  their  approach  [49].  This  discretization  procedure  can  be 
systematically  applied  to  many  problems,  and  is  named  the  discrete-convolutional 
method  of  moments  (DCMoM) .  The  DCMoM  is  actually  a  generalization  of  the 
approach  originally  used  by  Bojarski  [10],  and  is  briefly  described  in  Chapter 
5.  Ray  and  Mittra  [50]  and  Hurst  and  Mittra  [51]  used  a  DCMoM  discretization 
with  the  CGM  to  find  the  current  density  induced  on  large  plates  and  finite 
cylinders,  and  showed  that  the  procedure  could  be  efficiently  applied  to  treat 
matrix  equations  with  orders  in  excess  of  4000. 
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Van  den  Berg  used  the  CGM  and  several  related  algorithms  to  analyze  scat¬ 
tering  from  conducting  strips  and  dielectric  cylinders,  and  presented  a  con¬ 
vergence  rate  comparison  which  shows  that  these  algorithms  are  superior  to  a 
"simple"  iterative  algorithm  [52].  Bokhari  and  Balakrishnan  used  a  simple  itera¬ 
tive  algorithm  for  the  anaylsis  of  a  cylindrical  antenna,  and  reported  a  way  of 
always  ensuring  the  convergence  of  the  algorithm  by  the  use  of  an  appropriate 
amount  of  zero-padding  [53]. 

All  of  the  above  approaches  rely  on  the  presence  of  discrete-convolutional 
symmetries  in  the  matrix  equation  to  permit  a  reduction  in  the  necessary  com¬ 
puter  memory.  These  are  a  slightly  more  general  form  of  the  Toeplitz  and 
block-Toeplitz  symmetries  for  which  specialized  direct  methods  of  solution  have 
been  developed  [54l,  [55].  These  direct  algorithms  permit  the  same  storage 
reduction  as  iterative  techniques,  and  in  addition  may  be  more  efficient  than 
iteration.  However,  slightly  perturbed  Toeplitz  systems  and  Toeplitz  systems 
embedded  in  larger  matrix  equations  often  arise  in  practice,  and  are  not  easily 
treated  by  direct  algorithms.  For  instance,  the  pertinent  matrix  equations 
often  have  discrete-convolutional  symmetries  everywhere  except  along  the  main 
diagonal,  which  upsets  the  Toeplitz  form  of  the  system  without  changing  the 
significant  redundancy  features.  Examples  of  almost-Toepl itz  systems  are 
discussed  in  Chapters  6  and  7.  Historical ly ,  it  appears  that  the  primary  use  of 
iteration  in  electromagnetics  has  been  to  solve  equations  having  slightly  per¬ 
turbed  Toeplitz  symmetries. 

Because  of  the  reliance  on  discrete-convolutional  symmetries,  the  spe¬ 
cialized  techniques  discussed  above  are  only  optimum  for  a  limited  set  of 
problems.  These  include  geometries  which  can  be  described  by  surfaces  of 
constant  curvature  (such  as  flat  plates  or  right-circular  cylinders)  and 
geometries  which  lend  themselves  to  volume  discretization  (such  as  penetrable 
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dielectric  bodies).  Although  it  may  always  be  possible  to  treat  an  arbitrarily 
shaped  structure  with  a  volume  discretization,  this  approach  often  requires  a 
large  number  of  fictitious  unknowns  and  can  not  be  considered  optimum  for 
problems  other  than  penetrable  bodies.  In  other  words,  in  a  given  situation  it 
may  not  be  possible  to  cast  the  matrix  equation  into  discrete-convolutional 
form.  An  alternate  possibility  for  the  efficient  analysis  of  large  scatterers 
that  does  not  require  the  matrix  to  possess  discrete-convolutional  symmetries 
is  the  matrix  element  regeneration  (MER)  scheme  recently  proposed  by  Peterson, 
Sultan,  and  Mittra  [56]  -  [59].  Chapter  4  discusses  the  implementation  of  the 
MER.  Since  the  MER  does  not  rely  on  symmetries  in  the  system,  the  approach  is 
not  generally  as  efficient  as  the  specialized  techniques  discussed  above. 

The  above  discussion  has  emphasized  recent  trends  in  the  numerical  analysis 
of  electrically  large  scattering  problems..  The  purpose  of  this  report  is  to 
interpret  these  trends,  extend  the  state  of  knowledge  concerning  certain  aspects 
of  these  approaches  which  have  heretofore  gone  unreported,  and  effectively 
implement  similar  techniques  for  other  applications.  Specifically,  Chapter  2 
presents  several  iterative  algorithms  related  to  the  conjugate  gradient  method 
and  shows  how  these  algorithms  can  be  obtained  from  two  different  viewpoints. 
This  material  may  be  beneficial  to  anyone  attempting  to  modify  existing 
algorithms  to  improve  their  rate  of  convergence,  as  has  been  recently  proposed 
[60].  Chapter  3  illustrates  the  convergence  of  the  iterative  algorithms  of 
Chapter  2  when  applied  to  integral  equations  of  typical  electromagnetic  scat¬ 
tering  problems.  Chapter  4  begins  the  topic  of  efficient  implementation  of 
iterative  algorithms,  with  an  investigation  into  the  matrix  element  regeneration 
(MER)  approach.  Chapter  5  provides  an  in-depth  comparison  of  two  discretiza¬ 
tions  often  used  in  conjunction  with  iterative  methods,  the  discrete- 
convolutional  method  of  moments  (DCMoM)  and  the  spectral-domain  fast-Fourier 
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transform  (SDFFT)  procedures,  the  latter  denoting  the  discretization  used  with 
the  original  SIT  [11].  Chapters  6  and  7  concern  the  efficient  implementation  of 
iterative  approaches  for  the  problems  of  scattering  from  a  dielectric  body  and 
scattering  by  a  finite  hollow  conducting  cylinder.  Conclusions  and  suggestions 
for  future  work  are  given  in  Chapter  8. 

The  work  reported  within  is  intended  to  build  upon  the  previous  efforts  of 
others,  especially  R.  Kastner,  M.  Sultan,  and  S.  Ray.  Based  upon  the  foundation 
set  by  these  and  other  researchers,  the  application  of  iterative  methods  to 
numerical  electromagnetics  has  today  reached  a  stage  where  the  systematic  imple¬ 
mentation  of  a  variety  of  problems  is  we  1 1 -understood .  Appropriate  choices  may 
be  made  from  available  discretizations  and  available  iterative  algorithms,  in  a 
relatively  independent  manner.  As  improvements  are  made  in  either  of  these  two 
areas,  they  can  be  incorporated  into  the  existing  approaches.  This  study 
examines  several  aspects  of  iterative  computational  electromagnetics,  and  should 
be  considered  a  small  step  in  the  continuing  attempts  to  improve  the  efficiency 
of  numerical  techniques  for  solving  electromagnetic  scattering  problems. 
Preliminary  results  from  this  investigation  have  been  previously  published 
[56],  [57],  [59],  [61],  [62]. 
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2.  THE  CONJUGATE  GRADIENT  METHOD  AND  RELATED  ALGORITHMS 

2.1.  Introduction 

The  computer-aided  analysis  of  electrically  large  scatterers  is  often  based 
upon  an  iterative  solution  algorithm  [10]  -  [13].  In  this  chapter,  several 
iterative  "gradient"  methods  are  developed  and  presented  in  a  form  suitable  for 
this  type  of  analysis.  Of  central  importance  is  the  conjugate  gradient  method 
(CGM),  which  is  used  extensively  in  later  chapters.  In  the  literature,  the  CGM 
is  usually  presented  in  a  specialized  form,  which  is  useful  only  for  symmetric, 
positive  definite  systems  [30],  [31],  [33],  [63],  [64].  Here,  the  general 
algorithm  is  developed  from  two  different  perspectives.  In  addition,  several 
approaches  that  are  related  to  the  CGM  are  discussed. 

Initially,  the  CGM  is  developed  from  a  generalized  variational  procedure, 
i.e.,  the  minimization  of  an  error  functional.  The  CGM  can  also  be  developed 
from  an  orthogonal  expansion  process,  without  reference  to  a  functional. 

However,  it  will  be  apparent  that  these  two  constructs  are  actually  linked,  that 
is,  each  functional  has  associated  with  it  a  certain  type  of  orthogonal  expan¬ 
sion.  The  algorithms  reduce  to  the  process  of  generating  the  orthogonal  func¬ 
tions  and  finding  the  proper  coefficients  to  represent  the  desired  solution. 

It  suffices  to  consider  the  nonsingular  matrix  equation 

Lf  -  h  (2.1) 

where  L  denotes  a  linear  operator  (N  *  N  matrix),  f  is  an  unknown  function  to  be 
determined  (N  x  1  matrix),  and  h  is  a  given  right-hand  side  (also  an  N  x  1 
matrix).  An  inner  product  and  its  associated  norm  are  defined 


<f.g>  ”  gTf 


(2.2) 
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| |f | |  -  /<fTf>  (2.3) 

where  the  *T'  denotes  matrix  transpose  and  the  bar  complex  conjugation.  It  is 
desired  to  iteratively  generate  the  solution  f,  so  that  at  some  stage  of  the 
process  we  seek  an  estimate  of  the  solution  in  the  form 


f  =  f  ♦  a  p 

n  n— 1  nr r 


(2.4) 


where  fn_j  is  a  previous  estimate  of  the  solution,  pn  is  a  function  which  is  as 
yet  undetermined,  and  is  a  scalar  coefficient. 

2.2.  Functional  Minimization 

Assuming  a  function  pn  is  available,  the  coefficient  an  in  Equation  (2.4) 

can  be  chosen  so  that  f^  is  a  "better"  estimate  of  the  solution  than  ffj_1 ,  in 

the  sense  that  f  minimizes  some  measure  of  the  error.  For  instance,  the 
n 

quadratic  functionals 


F.(f  )  =  | | Lf  -  h| 
In  1  1  n  1 


F,(f  )  -  I  If  -  f 

2  n  11  n 


(2.5) 

(2.6) 


provide  two  possible  definitions  of  the  error  in  f  .  The  functional 

n 


F,(f  )  =  <Lf  -  h,  f  -  f  > 
J  n  n  n 


(2.7) 


can  also  be  used  to  define  the  error  in  f^,  although  a  useful  iterative  proce¬ 
dure  can  be  based  on  only  if  L  is  self-adjoint  and  positive  definite  in  the 

sense  of  Griffel  (65).  Each  of  the  above  functionals  gives  rise  to  a  different 
iterative  procedure  for  the  solution  of  Equation  (2.1);  other  functionals  could 
also  be  used  and  would  lead  to  other  algorithms. 
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Functional  F^,  as  defined  in  Equation  (2.5),  ia  associated  with  the  stan¬ 
dard  form  of  the  CGM  for  arbitrary  linear  operators.  For  a  given  pn>  the  coef¬ 
ficient  of  Equation  (2.4)  that  minimizes  Fj  is  found  by  setting  the  first 
variation  of  F^  with  respect  to  to  zero,  which  is  equivalent  to 

2Re{5(T")(<Lf  ,  -  h,Lp  >  +  a  I  |Lp  |j2)}  =  0  (2.8) 

1  a  n- 1  n  n 1  1  rn M  1 

for  arbitrary  variation  5(1P) .  This  is  satisfied  for 


a 

n 


I tLpn | (2 


(2.9) 


where  for  convenience  we  define 


R  =  Lf  -  h  (2.10) 

n  n 

If  it  is  desired  to  simultaneously  find  two  coefficients  which  minimize  F^f^) 
where 


the  above  process  can  be  extended  to  produce 

5  _  '^n-l'^HLqJI2  +  <Rn-l>LV<LVLPn> 

llLPnll2  I  |Lqn  II2-  l<LPn>L£in>i2 


(2.11) 


(2.12) 


7  .  -<Rn-PLVNLPnM2  +  <Rn-l,LPn><LPn,LV  (2  ^ 

ll^Pnl!2  IlLqJI2  -  |<Lpn.Lqn>|2 

From  Equations  (2.12)  and  (2.13)  it  is  apparent  that  if  the  functions  p^  and  qn 
satisfy 
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<Lpn,Lqn>  *  0 


(2.14) 


the  expressions  for  the  coefficients  decouple  and  reduce  to  the  form  presented 
in  Equation  (2.9).  This  indicates  that  if  a  sequence  of  functions  satisfying 
Equation  (2.14)  is  used  to  expand  the  solution  f  of  Equation  (2.1),  the  coef¬ 
ficient  of  each  function  may  be  found  independently.  Under  these  conditions, 
the  process  cannot  be  improved  upon  by  simultaneously  computing  several  coef¬ 
ficients  . 

It  is  also  of  interest  to  consider  a  solution  estimate  of  the  form 


f  =  f  ,  +  a  ( p  +■  8  q  ) 
n  n-1  n  Fn  nMn 


(2.15) 


In  the  manner  outlined  above,  the  8  that  minimizes  F,(f  )  assuming  optimum  a 

n  in  r  n 

is  given  by 


s  I  lLP„i  I  - 

"  <R,.1.‘-pn>llMtlll2  -  Vl'W'V 


(2.16) 


Suppose  that  p^  and  qn  are  arbitrary  functions  which  do  not  satisfy  Equation 
(2.14).  Suppose  also  that  q^  has  already  been  used  in  the  expansion  process,  so 


f  .  =  f  „  +  a  .  q 
n-1  n-2  n-l^n 


(2.17) 


where  is  found  according  to  Equation  (2.9).  It  immediately  follows  that 


<R  , ,Lq  >  *  <R__,  +  a  . Lq  ,Lq  >  *  0 
n  I  n  n-£  n-i  n  ’n 


(2.18) 


_<LPn,Lqn> 


(2.19) 
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and 


<Lqn,L(pn  ♦  0nqn)>  =  0  (2.20) 

Thus,  the  process  of  selecting  coefficients  that  minimize  the  functional  F^  is 
optimized  when  functions  satisfying  Equation  (2.14)  are  used.  If  arbitrary 
functions  are  used,  and  their  coefficients  adjusted  simultaneous ly  to  produce  a 
minimum  F^,  the  process  will  adjust  these  coefficients  in  order  to  produce  func¬ 
tions  satisfying  Equation  (2.14). 

2.3.  The  Conjugate  Direction  Method 

The  previous  section  presented  a  procedure  for  finding  coefficients  of 
trial  functions  that  minimize  an  error  functional,  for  the  solution  of 
Equation  (2.1).  It  was  seen  that  a  sequence  of  functions  (pnf  satisfying 

<Lpi,Lp->  =0  i  *  j  (2.21) 

arose  naturally  in  connection  with  the  error  functional  Fj  defined  in 
Equation  (2.5).  If  we  assume  the  existence  of  such  a  sequence,  and  note  that  in 
a  finite  dimensional  vector  space  N  of  these  p-functions  will  span  the  space, 
the  solution  can  be  written 

f  =  fQ  +  alPl  +  +  Vn  +  "•  +  VN  '  (2-22) 

where  fg  is  an  arbitrary  initial  estimate  of  the  solution.  Because  of  the 
orthogonality  expressed  in  Equation  (2.21),  each  coefficient  can  be  found  inde¬ 
pendently  as 


~<R0,Lpn> 
I  I Lpn  I  1 2 


a  a 
n 


(2.23) 
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where  is  defined  in  Equation  (2.10).  Note  that  Equation  (2.23)  and  Equation 
(2.9)  appear  to  differ;  it  will  be  shown  below  that  they  are,  in  fact,  equal. 
From  the  above  relationships,  it  is  apparent  that 


R  3  8.  a.L-  +  ...  +  a  Lp 
n  0  1  r  1  n  rn 

and  recursive  relationships  are  given  as 


(2.24) 


R  **  R  .  +  a  Lp 
n  n-1  n  *n 

f  =  f  ,  +  a  p 
n  n-1  nrn 

The  "residual  norm"  defined  by 

_  I  I Rn  1  I  _  l|Lgn  -  hi  | 
=  INI  ‘  f|h|| 


(2.25) 

(2.26) 


(2.27) 


is  one  indication  of  the  accuracy  of  the  solution  at  the  n-th  step  in  the  expan¬ 
sion  process.  Because  of  the  orthogonality  of  the  p-functions,  Equation  (2.24) 
can  be  used  to  show  that 


and,  thus,  the  residual  norm  decreases  monotonical ly  as  n  increases.  Of  course, 
the  residual  norm  is  directly  proportional  to  the  error  functional  F^  discussed 
above,  and  must  decrease  at  each  step  if  a  minimization  is  performed.  From 
Equations  (2.21),  (2.23),  and  (2.24)  we  can  readily  deduce  that 

<R0»Lpm>  n  <  m 

<Rn,LP.>-<  (2'29) 

o  n  2.  m 

V 
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This  relationship  shows  the  equivalence  of  Equations  (2.9)  and  (2.23)  and  will 
prove  useful  below. 

Functions  that  satisfy  Equation  (2.21)  are  said  to  be  mutually  conjugate 
with  respect  to  the  operator  L^L,  where  LA  is  the  adjoint  to  L  defined  by 

<LAR,P>  =  <R,LP>  (2.30) 


The  above  process  of  expanding  a  solution  in  conjugate  functions  is  known  as  the 
"conjugate  direction  method,"  after  Hestenes  and  Stiefel  [63].  The  conjugate 
direction  method  does  not  specify  a  procedure  for  the  generation  of  mutually 
conjugate  functions.  We  now  consider  such  a  procedure. 

2.4.  The  Conjugate  Gradient  Method  (CGM) 

An  approach  that  is  based  on  the  conjugate  direction  method  and  includes 
an  algorithm  for  the  generation  of  mutually  conjugate  p-functions  was  developed 
by  Hestenes  and  Stiefel  and  named  the  "conjugate  gradient”  method  [63].  The 
process  begins  with  the  choice 

Pj  =  -lVq  (2.31) 

Subsequent  functions  are  found  from 


(2.32) 


where  the  scalar  coefficient  8  is  chosen  to  ensure 

n 


<1*Al'pn’pn+l>  ’  0 


(2.33) 


It  will  be  shown  that  due  to  the  special  nature  of  the  R-functions,  this  process 
yields  a  mutually  conjugate  set.  First,  some  preliminary  relationships  are 
presented . 
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Based  upon  Equation  (2.32),  we  write 

<o  L^R  >  =  -<LAR  LAR  >  +  $  <p  ,L^R  > 
^Pn+1  *L  vlj  n>L  m  n  *n  ’  m 


(2.34) 


From  Equation  (2.29),  the  first  and  last  inner  products  in  Equation  (2.34) 
vanish  for  m  >  n,  leaving 


<LAR  ,LAR  >  =  0  m  *  n 


(2.35) 


Equations  (2.32)  and  (2.29)  can  be  combined  to  yield 


<lar  ,p  aI>  =  -<lar  ,lar  > 

n,Kn+l  n’  n 


(2.36) 


It  follows  from  Equation  (2.23)  that 

llL\-ll|2 

a  =  - s — 

I|M>„II 

From  Equation  (2.25),  we  have 


(2.37) 


=  LAR  .  +  a  lAp 
n  n-1  n  n 


(2.38) 


An  inner  product  between  L  R^  and  Equation  (2.38)  leads  to  the  result 


r 


llL\ 


*\ 


-NlArJI2 


n  m 


m  =  n  -  1 


otherwise 


(2.39) 


Using  Equation  (2.39)  with  ra  =  n,  we  find  the  value  of  £?n  from  Equations  (2.32) 


and  (2.33)  to  be 
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lL\l 


|lAr  ,|j2 

i  n-l  11 


(2.40) 


To  see  that  Bn  guarantees  the  proper  orthogonality  between  functions  when 
Equation  (2.21)  is  not  enforced  explicitly,  consider  the  above  process  as  an 
iterative  method.  Given  some  "initial  guess"  fQ,  we  generate 

Rq,  pj ,  Rj ,  P2 ,  ...  in  the  prescribed  manner.  Since  Equation  (2.33)  is  expli¬ 
citly  enforced  we  know  that 


<larq,lAr2>  =  0 
<LARrLAR2>  =  0 


What  remains  is  the  validity  of 


<LAbp1  ,p3>  =»  0 


By  direct  computation  p^  can  be  written  as 


-IIlV 


L  R. 


i-o  I  |lar,  1 12 


(2.41) 

(2.42) 

(2.43) 


(2.44) 


(2.45) 


and  it  follows  that 


.  .  -  2  <LALp, ,LAR.> 

<Ap1,P3>.-|M;II 


(2.46) 

But,  by  the  relationship  established  in  Equation  (2.39),  the  above  reduces  to 


o 


(2.47) 
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Thus,  in  an  inductive  fashion  we  see  that  the  functions  generated  by  the  above 
procedure  do  in  fact  satisfy  the  assumed  orthogonality  properties  of  the  con¬ 
jugate  direction  method. 

The  conjugate  gradient  method  is  given  in  its  entirety  as  follows: 


Initial  Steps 


Guess  f 

o 


Ro  h 

p,  -  -i\ 


Iterate  n  -  1,2,... 


f 

n 


"R\-,n2 

iKii2 

f  ,  ♦  Cl  p 

n-1  rvn 


<Vl’LP»> 

II1P.II2 


R  *  Lf  -  h  *  R  ,  +  a  Lp 
n  n  n-i  n  rn 

1 1 12 

6  — 2_ 

n  i  ii\.i  1 1 

p  xl  =  -LAR  +  s  p 
*n+l  n  n^n 

A  discussion  of  the  theoretical  convergence  rate  of  the  CGM  and  the  con¬ 
vergence  in  practice  is  presented  in  Chapter  3.  The  remainder  of  this  chapter 
considers  alternate  approaches  that  are  related  to  the  algorithms  discussed 


above , 
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Ynk  =  0  k  <  n  -  l  (2.52) 

Because  of  this,  the  CGM  generates  “  orthogonal  functions  recursively, 
without  the  need  to  permanently  store  the  p-functions  and  R-functions  in  com¬ 
puter  memory.  This  process  is  summarized  as  follows,  based  on  a  given  function 


(2.53) 
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R 

n 


P 


n+1 


Vi  + 


l^„-r2  Lpn 


I  |LPr 


-lV  .  -ifol1.,  , 


”  nA-iii 


2  n 


(2.54) 


(2.55) 


As  an  illustration  of  a  scheme  that  produces  the  same  R-functions  but  nor¬ 
malizes  the  p-functions  differently,  consider 


R 

n 


pn+l 


LP„ 


Vi  + 


lLPnl 


LARn 

TiSj? 


(2.56) 


(2.57) 


(2.58) 


Other  normal izations  could  also  be  used,  and  would  be  equivalent  to  the  CGM 
procedure . 

2.6.  Minimization  in  the  Domain  of  the  Operator  (PGM) 

The  CGM  as  presented  above  was  shown  to  minimize  the  functional  defined 
in  Equation  (2.5)  at  each  step  of  the  procedure.  Functional  F^  measures  the 
error  in  a  given  estimate  of  the  solution  as  seen  in  the  range  of  the  operator 
L.  CGM-like  algorithms  can  be  developed  for  many  other  functionals  [66],  [67]. 
For  instance,  a  similar  algorithm  can  be  based  upon  functional  F2  defined  in 
Equation  (2.6),  and  repeated  for  convenience  below: 

W  -  Mf  -  fjl2  (2.59) 


F2  provides  an  error  measure  made  in  the  domain  of  the  operator  L.  If  a 
solution  is  sought  in  the  form 
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f  =  f  .  +  a  q 
n  n-1  nMn 


(2.60) 


the  value  of  the  coefficient  an  which  minimizes  F2  is  given  by 

<f  -  f  ,  ,q  > 
n-1  ,Mn 


a  = 
n 


(2.61) 


Since  we  seek  a  solution  to  Equation  (2.1),  we  do  not  know  f  at  this  point. 
However,  we  can  use  the  identity 


(l*1)a  =  (lV1 


(2.62) 


and  the  previous  definition 


R  =  L(  f  -  f) 
n  n 


(2.63) 


to  rewrite  Equation  (2.61)  as 


(2.64) 


which  is  shown  below  to  be  a  useful  form. 

At  this  point,  consider  a  solution  to  E- urtion  (2.1)  in  the  form 


f  =  fo  +  Vi  +  •••  +  Vn 


(2.6b) 


where  the  functions  (qn)  span  the  N-dimensional  space  and  satisfy  the  orthog¬ 
onality  requirement 


<qL.qj>  ”0  i  *  j 


(2.66) 


It  is  apparent  that  the  coefficients  are  given  by 

-<vaA,'V 


2 


(2.67) 
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without  regard  to  minimizing  any  specific  functional.  By  analogy  with  the  CGM 
development  described  previously,  we  conclude  that  the  functional  is  asso¬ 
ciated  with  an  orthogonal  expansion  of  the  type  appearing  in  Equation  (2.66). 
Note  that,  at  this  point  in  the  development,  (LA)  1  is  not  known.  However,  by 


making  the  choice 
<11  =  ~L\ 

and 

q  ,  *  — LAH  +  0  q 

Hn+1  n  nHn 

we  can  construct  an  algorithm  that  does 
A  ~  1  - 1 

(L  )  or  L  .  Furthermore,  it  can  be 

R  =  R  .  +  a  Lq 
n  n-1  n 

and 


(2.68) 

(2.69) 

not  require  the  explicit  use  of  either 
shown  that 

(2.70) 


<R  ,R  >  =  0 
n  m 


n  *  cb 


(2.71) 


n- 1 


(2.72) 


8. 


(2.73) 


n-  1 


by  carrying  out  a  development  which  parallels  that  given  for  the  CGM  case  in 
Section  2,4.  We  note  that,  strictly  speaking,  the  minimization-in-the-domain 
algorithm  is  not  a  conjugate  direction  method,  because  the  orthogonality  between 
expansion  functions  is  simple  orthogonality,  as  expressed  in  Equation  (2.66). 
Perhaps  it  is  clearer  to  name  this  method  an  "orthogonal  gradient"  method  (OGM). 
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The  complete  OGM  algorithm  is  summarized  below: 


Initial  Steps 
Guess  f 


o 


R 

q 


0 

1 


Lf  - 


o 


h 


Iterate  n  =  1,2,... 


a 

n 


n-1 


2 


f  =  f  ,  +  a  q 

n  n- 1  n4n 

Rn  =  Lfn  "  h  =  Rn  !  +  a„L<l„ 

n  n  n-I  n  ^n 


q  .  =  -L^R  +  8  q 
n+1  n  n^n 

2.7.  A  Gradient  Algorithm  Using  Approximate  Inverse  Operators  (AIGM) 

From  the  above  discussion  of  the  CGM  and  the  OGM,  it  is  clear  that  the 
algorithms  generate  a  sequence  of  expansion  functions  and  then  find  their  coef¬ 
ficients  to  minimize  some  measure  of  the  error.  Furthermore,  the  CGM  and  OGM 
both  generate  expansion  functions  which  are  mutually  orthogonal,  and  thus, 
theoretically  terminate  in  a  finite  number  of  iterations.  These  algorithms  are 
in  some  sense  equivalent  to  generating  an  inverse  operator  to  a  prescribed 
degree  of  accuracy.  In  many  applications,  an  approximate  inverse  operator  is 
readily  available,  but  not  accurate  enough  to  produce  a  meaningful  solution 
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directly.  In  this  section,  we  consider  a  gradient  algorithm  that  incorporates 
an  approximate  inverse  operator  in  order  to  systematically  improve  an  estimate 
of  the  solution.  The  rate  of  convergence  of  this  type  of  algorithm  depends  on 
the  accuracy  of  the  approximate  inverse  operator;  unfortunately,  there  is  no 
theoretical  guarantee  that  convergence  will  occur  in  a  finite  number  of  itera¬ 
tions.  This  procedure  is  a  generalization  of  an  approach  used  by  van  den  Berg 
[52]. 

If  the  standard  form  of  the  CGM  is  modified  so  that 


(2.74) 

(2.75) 


where  L 


is  the  approximate  inverse  operator,  and  if  is  chosen  so  that 


<Lpn+1,Lpn>  *  0  (2.76) 

the  resulting  algorithm  provides  a  systematic  way  to  generate  the  solution  to 
Equation  (2.1).  The  motivation  for  this  choice  stems  from  the  fact  that 


which  follows  from  Equation  (2.10).  Thus,  if  L  ^  is  the  exact  inverse,  the  pro¬ 
cess  terminates  after  the  generation  of  p^  according  to  Equation  (2.74).  Note 
that  the  above  choice  of  8n  does  not  ensure  that  the  sequence  (pn)  is  mutually 
-  orthogonal,  as  it  did  with  the  CGM.  The  orthogonality  relationships  that 
were  described  for  the  CGM  do  not  apply  to  this  algorithm,  and  a  finite  number 
of  the  p-functions  will  not,  in  general,  span  the  solution  space.  Even  if  they 
did,  their  coefficients  are  not  adjusted  accordingly,  but  rather  with  the  form¬ 
ula  of  Equation  (2.9).  We  emphasize  that  the  usefulness  of  this  approach, 
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which  we  name  the  "approximate  inverse  gradient"  method  (AIGM),  is  directly  pro¬ 
portional  to  the  accuracy  of  L  ^ 

The  complete  AIGM  algorithm  is  given  as  follows: 


Initial  Steps 
Guess  f 


-  h 


R 


0 


Iterate  n  =  1,2,... 


IlLpJI2 


E„  •  f„-i  *  V. 


R  -  Lf  -  h  =  R  ,  +  aLp 
a  n  n-L  n  rn 


<LL_IR  ,Lp  > 
n  *n 


lLPn 


p  ,  *  -L  ^  R  8  p 
*n+l  n  n*n 


This  AIGM  is  based  on  minimizing  functional  defined  in  Equation  (2.5). 
Similar  algorithms  could  be  based  on  other  functionals. 

2.8.  Summary 

This  chapter  has  presented  several  procedures  that  can  be  used  for  the 
iterative  solution  of  matrix  equations,  with  the  intent  of  emphasizing  the 
important  aspects  of  their  development.  The  CGM  algorithm  is  the  common  form 
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of  the  conjugate  gradient  method  for  arbitrary  matrix  equations,  and  is  some¬ 
times  known  as  the  minimizat ion-in-the-range  method.  The  OGM  is  a  different 
algorithm  which  is  sometimes  known  as  the  minimization-in-the-domain  version  of 
the  CGM.  The  AIGM  is  a  procedure  that  is  based  on  a  minimizat ion-in-the-range 
procedure,  and  allows  the  systematic  use  of  an  approximate  inverse  operator.  An 
important  feature  of  these  three  algorithms  is  that  they  will  not  diverge,  as 
simple  iterative  methods  such  as  the  Jacobi  and  Gauss-Seidel  algorithms  [30]  - 
[34]  occasionally  do  when  applied  to  problems  of  practical  interest.  A 
discussion  of  the  performance  of  the  CGM,  OGM,  and  AIGM  is  presented  in  Chapter 
3,  along  with  some  theoretical  results  concerning  the  convergence  of  the  CGM. 
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3.  CONVERGENCE  OF  ITERATIVE  ALGORITHMS  WHEN  USED  FOR 
ELECTROMAGNETIC  SCATTERING  APPLICATIONS 

3.1.  Introduction 

The  previous  chapter  presented  several  iterative  algorithms  but  omitted  a 
discussion  of  their  convergence.  This  chapter  provides  an  illustration  of  the 
convergence  for  several  types  of  problems  arising  in  electromagnetic  scattering. 
A  variety  of  examples,  intending  to  represent  typical  performance,  are  pre¬ 
sented.  While  the  conclusions  from  this  study  are  not  necessarily  applicable  to 
every  problem  that  may  arise,  they  will  serve  as  useful  guidelines  in  the 
characterization  of  the  iterative  algorithms  for  many  practical  electromagnetic 
scattering  problems.  Furthermore,  it  appears  that  a  familiarity  with  the  per¬ 
formance  of  the  iterative  algorithms  is  useful  in  a  different  way:  it  is  helpful 
in  identifying  situations  where  the  numerical  modeling  is  inadequate  to  produce 
a  meaningful  solution. 

Chapter  2  presented  three  algorithms,  the  conjugate  gradient  method  (CGM), 
the  orthogonal  gradient  method  (OGM)  and  the  approximate-inverse  gradient  method 
(AIGM) .  Here,  implementations  of  all  three  of  the  algorithms  are  considered,  the 
majority  being  examples  of  the  CGM.  The  rate  of  convergence  of  the  OGM  is  very 
similar  to  that  of  the  CGM,  and  in  general,  the  performance  of  the  AIGM  depends 
on  the  particular  approximate-inverse  operator  in  use.  In  all  cases  discussed 
here,  the  approximate  inverse  is  obtained  by  inverting  the  main  diagonal  of  the 
matrix.  This  specific  implementation  of  the  AIGM  appears  to  converge  slower  than 
the  CGM  for  most  of  the  examples  presented  herein,  although  it  is  a  more  general 
implementation  of  the  AIGM  than  that  used  by  van  den  Berg  [52].  Van  den  Berg 
presented  an  example  for  which  a  more  accurate  approximation  to  the  inverse  was 
available,  enabling  the  AIGM  to  consistently  outperform  the  CGM  [52]. 
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Primarily,  we  are  concerned  with  the  solution  of  the  matrix  equation 
arising  from  a  moment-method  discretization  of  an  integral  (or  integro- 
differential)  equation,  as  discussed  by  Harrington  [25]  and  Wilton  and  Butler 
[68].  The  examples  used  for  illustration  are  all  equations  representing  electro¬ 
magnetic  scattering  problems.  The  formulation  used  for  each  example  is  not 
discussed  in  detail  here,  but  is  available  in  later  chapters  or  elsewhere  in  the 
literature.  Furthermore,  at  this  point  we  are  not  concerned  with  the  implemen¬ 
tation  of  these  algorithms  for  the  efficient  treatment  of  large-order  systems  of 
equations.  This  topic  is  reserved  for  later  chapters.  Thus,  most  of  the  examples 
presented  here  are  based  upon  matrix  equations  of  relatively  small  order,  i.e., 
below  order  150. 

The  convergence  behavior  of  the  iterative  algorithms  will  be  illustrated 
via  the  residual  norm 


Lfk  -  h|  I 

Tmi 


(3.1) 


Experimentation  with  this  quantity  has  shown  it  to  be  reliable  provided  that  the 
matrix  L  is  not  badly  conditioned.  Normally,  a  criterion  of 


\  <  10 


(3.2) 


is  adequate  to  ensure  several  digits  of  accuracy  in  the  solution.  The  condition 
of  the  matrix  L  has  a  significant  effect  on  the  performance  of  the  iterative 
algorithms,  as  will  be  discussed  in  Section  3.2. 

All  of  the  iterative  algorithms  require  the  user  to  provide  an  initial 
estimate  or  "guess"  of  the  solution.  In  all  cases  examined  here,  a  zero  estimate 
is  used.  All  of  the  algorithms  considered  here  will  converge  in  theory  for 
arbitrary  initial  estimates,  a  feature  not  shared  by  many  of  the  iterative 
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algorithms  used  in  the  past.  In  addition,  the  rate  of  convergence  is  not  a 
strong  function  of  the  initial  guess,  as  is  indicated  by  the  theoretical  rates 
of  convergence  presented  in  Section  3.2. 

A  desirable  feature  of  any  algorithm  for  the  solution  of  electromagnetic 
scattering  problems  is  the  ability  to  solve  Equation  (2.1)  for  many  different 
right-hand  sides.  Specifically,  it  may  be  necessary  to  generate  solutions 
corresponding  to  many  different  orientations  of  the  scatterer  and  excitation  in 
space,  all  of  which  can  be  described  by  the  same  matrix  operator  L.  Some  com¬ 
putational  savings  could  arise  by  simultaneously  generating  several  solutions 
with  the  same  set  of  expansion  functions  (the  p-functions  discussed  in  Section 
2.3).  This  approach  is  examined  in  Section  3.7,  but  it  does  not  appear  to  yield 
the  expected  savings  due  to  properties  of  the  iterative  algorithms  in  use. 

3.2.  Aspects  of  the  Theoretical  Convergence  of  the  CGM 

From  the  discussion  of  the  CGM  in  Chapter  2,  it  is  apparent  that  in  the 
absence  of  any  machine  errors  in  the  various  ca leu lat ions ,  the  CGM  will  ter¬ 
minate  in  N  steps  or  less  for  an  order  N  system.  In  fact,  the  number  of  itera¬ 
tions  required  for  the  solution  of  Equation  (2.1)  is  normally  equal  to  the 
number  of  independent  eigenvalues  in  the  matrix  iA  [34],  [69],  Furthermore,  if 

f  is  the  estimate  of  the  solution  f  produced  by  the  CGM  after  n  iterations,  f 
n  J  n 

must  satisfy 

Hfn  ‘  fll  <  Mfn  ‘  fM  .  "  >  m  (3.3) 


Thus,  the  estimates  produced  by  the  CGM  converge  raonotonical ly  to  the  solution. 
Equation  (3.3)  may  be  obtained  by  considering  the  general  form  of  Equation 


(2.45) 
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P 


n 


n-1  LAR£ 

i-0  IlL^.H2 


From  Equation  (3.4)  and  the  orthogonality  of  the  set 
Equation  (2.35),  we  have 


(3.4) 


(LAR.  }  as  expressed  in 


<Pi  .Pj>  >.  0 

From  the  definition  of  f^,  it  follows  that 


(3.5) 


f 

n 


f 

m 


n 


I 

i=m+l 


a.  p . 


n  >  m 


(3.6) 


Note  that  the  {a^}  are  nonnegative  by  Equation  (2.37).  Equations  (3.5)  and  (3.6) 
can  be  combined  to  yield  the  inequality 


<fn  -  -  fn>  >  0  N  >  n  >  m 

n  m  N  n  — 


Finally,  Equation  (3.7)  can  be  combined  with  the  identity 


(3.7) 


Hf»  *  f  M2  *  Mfm  -  «nM2  *  Hfn  -  f  I  I2  -  2Re{<fn  -  fn,f  -  fn>}  (3.8) 

to  prove  Equation  (3.3).  It  is  interesting  to  note  that  both  the  CGM  and  the  OGM 
algorithms  produce  solution  estimates  which  satisfy  Equation  (3.3);  the  AIGM 
algorithm  does  not.  Of  course,  the  presence  of  round-off  errors  may  invalidate 
Equation  (3.3)  to  some  extent. 

Because  round-off  errors  degrade  the  finite-step  termination  property  of 
the  CGM,  in  practice  it  is  often  considered  a  purely  iterative  algorithm.  In 
this  context,  an  upper  bound  on  its  rate  of  convergence  is  given  as  [67),  [69] 


33 


lfn  - 

lfo  -  fl 


r<  -  i 


<  2  — - 

/<  +  1 


(3.9) 


where  <  is  the  condition  number  of  iA.  (the  ratio  of  the  maximum  to  minimum 
eigenvalues  of  the  matrix  iA) .  In  addition,  it  has  been  observed  that  approxima¬ 
tely  n/ic  iterations  are  normally  required  to  ensure  accuracy  to  n  decimal  pla¬ 


ces  [69].  Data  to  follow  support  the  conclusion  that  the  dependence  of 


convergence  rate  on  the  number  of  decimal  places  is  approximately  linear. 

These  approximate  rates  of  convergence  are  independent  of  the  right-hand 
side  of  the  system  and  the  initial  estimate  of  the  unknown,  traits  which  were 
usually  observed  in  the  numerical  examples  examined  throughout  this  study.  In  an 
effort  to  reduce  the  amount  of  iterative  computation  as  much  as  possible,  appre¬ 
ciable  effort  is  often  expended  to  produce  a  "good"  initial  estimate  of  the 
solution.  For  instance,  asymptotic  solutions  to  electromagnetic  scattering 
problems  are  often  used  as  initial  estimates  in  iterative  solutions  [44).  Based 


upon  the  above  observations,  an  initial  estimate  will  probably  not  affect  the 
rate  of  convergence,  only  the  distance  from  the  starting  point  to  the  actual 
solution.  Since  the  rate  of  convergence  is  roughly  linear  with  the  required 
number  of  decimal  places  of  accuracy,  an  initial  estimate  can  be  considered 
"good"  only  if  it  reduces  to  the  initial  value  of  the  residual  norm  by  a  signi¬ 
ficant  amount  below  unity.  (Unity  is  the  value  of  SIQ  produced  by  an  initial  guess 
of  zero.)  As  an  example,  suppose  that  it  is  desired  to  produce  a  solution  with 

residual  norm  V  =  10  4.  An  initial  estimate  of  the  solution  which  produces 
-2 

=>  10  should  reduce  the  amount  of  iterative  computation  to  half  of  what  would 
have  been  required  if  an  initial  estimate  of  zero  had  been  used.  Similarly,  an 
initial  estimate  with  .V  *  10  ^  should  reduce  the  required  iterative  computation 
by  one  quarter,  and  so  on. 
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Whether  considered  a  finite-step  method  or  a  purely  iterative  method,  the 
rate  of  convergence  of  the  CGM  depends  on  the  eigenvalues  of  the  matrix.  The  NxN 
matrix  produced  by  a  moment-method  discretization  of  an  integral  equation  can  be 
considered  an  approximation  to  the  integral  operator.  Properties  of  the  integral 
operator,  such  as  its  eigenvalue  spectrum,  are  projected  onto  the  matrix.  In 
general,  both  the  accuracy  of  the  solution  to  the  system  and  the  accuracy  of  the 
eigenvalue  representat ion  will  be  linked  to  the  order  of  the  discretization,  and 
should  improve  as  the  discretization  is  refined. 

Suppose  that  a  given  integral  equation  is  discretized  via  the  moment  method 
to  yield  an  N-th  order  matrix  equation,  which  is  then  solved  by  the  CGM.  Let  us 
also  suppose  for  the  moment  that  the  discrete  system  we  obtain  will  "converge" 
in  some  sense  to  the  integral  equation  as  the  order  of  the  discretization  is 
increased.  If  a  solution  to  the  N-th  order  system  is  found  in  n  steps,  where  n 
is  much  smaller  than  the  order  of  the  system,  the  theoretical  property  of  the 
CGM  suggests  that  there  are  n  dominant  eigenvalues  in  the  spectrum  of  the  asso¬ 
ciated  integral  operator.  A  further  check  on  this  result  can  be  obtained  by  re¬ 
discretizing  the  integral  equation  to  yield  a  larger-order  matrix  equation,  for 
which  the  CGM  should  again  converge  in  n  steps.  This  behavior  is  typical  for 
problems  where  the  matrix  is  able  to  represent  the  eigenvalue  structure  of  the 
integral  operator,  as  illustrated  in  the  examples  to  follow. 

Suppose,  on  the  other  hand,  convergence  is  slow  in  that  n=N.  If  the  CGM  is 
then  applied  to  a  larger  system  arising  from  a  refinement  of  the  discretization, 
and  many  additional  iterations  are  necessary  for  convergence,  it  is  a  clear 
indication  that  the  original  matrix  is  not  able  to  represent  the  eigenvalue 
structure  of  the  integral  operator.  Under  such  conditions  the  particular  matrix 
equation  is  probably  not  an  accurate  model  of  the  integral  equation,  in  that  the 
level  of  discretization  is  inadequate  to  represent  the  problem.  Alternatively, 
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it  is  possible  that  the  formulation  on  which  the  matrix  equation  is  based  may  be 
invalid,  which  could  cause  a  significant  change  in  the  typical  behavior  of  the 
iterative  algorithm.  An  example  of  the  latter  situation  arises  when  an  integral 
equation  exhibits  non-unique  solutions,  such  as  occur  in  connection  with  the 
"internal  resonance"  problems  discussed  briefly  in  Section  3.3. 

Although  the  preceding  remarks  are  somewhat  speculative,  numerical  results 
support  the  notion  that  the  CGM  provides  feedback  which  can  be  useful  for  iden¬ 
tifying  situations  where  the  accuracy  of  the  numerical  modeling  process  appears 
questionable.  This  behavior  is  illustrated  throughout  the  remainder  of  this 
chapter.  However,  not  every  example  will  conform  to  the  same  type  of  behavior. 
User- fami 1 iar ity  with  the  typical  behavior  of  the  algorithm  for  the  particular 
type  of  problem  involved  is  essential  before  unusual  behavior  can  be  identified; 
even  then  the  above  ideas  are  not  absolutes  and  the  behavior  may  be  misleading. 
If  the  operator  is  ill-conditioned,  severe  round-off  errors  will  disrupt  the 
finite  step  termination  property  of  the  CGM  and  mask  the  convergence  behavior. 
Normally,  under  such  conditions  the  slow  convergence  of  the  algorithm  is  a  clear 
indication  that  the  system  is  poorly  conditioned,  which  itself  may  be  useful 
informat ion . 

Round-off  errors  occur  to  some  degree  whenever  machine  computations  are 
performed,  and  thus  mandate  the  choice  of  algorithms  which  minimize  their 
effects  whenever  possible.  Round-off  errors  affect  the  CGM,  which  generates 
expansion  functions  recursively  to  satisfy  the  orthogonality  of  Equation  (2.21), 
by  degrading  the  actual  orthogonality  of  these  functions.  If  the  round-off 
errors  are  severe,  the  functions  may  not  even  approximately  span  the  N- 
dimensional  solution  space.  Due  to  this  and  the  additional  fact  that  the  coef¬ 
ficients  of  the  expansion  functions  are  also  incorrect  due  to  round-off  errors, 
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the  algorithm  may  not  converge  in  practice.  Simon  illustrates  the  loss  of  orthog¬ 
onality  in  recursively  generated  functions  for  a  similar  algorithm  [70].  It  has 
been  noted  that  the  stability  of  the  CGM  to  round-off  errors  can  be  very  low 
[64]. 

Although  it  is  possible  to  examine  the  theoretical  stability  of  numerical 
algorithms  to  round-off  errors,  it  is  much  easier  and  perhaps  more  meaningful  to 
experiment  with  different  procedures  and  observe  their  performance  in  practice. 
Round-off  errors  did  not  appear  to  present  a  problem  for  the  examples  taken  from 
electromagnetic  scattering  problems  illustrated  in  the  remainder  of  this  chapter. 
These  examples  were  generated  with  the  aid  of  a  CDC  Cyber-175  computer,  which 
uses  a  single-precision  word  length  of  60  bits.  However,  when  poorly  conditioned 
systems  were  constructed  for  the  sole  purpose  of  testing  the  CGM,  the  algorithm 
was  not  able  to  solve  them  to  any  reasonable  degree  of  accuracy.  This  behavior 
had  prompted  past  research  into  CGM-like  algorithms  with  better  numerical  stabil¬ 
ity,  and  several  of  these  are  discussed  by  Stoer  [69].  A  different  avenue  of 
approach  has  been  the  use  of  preconditioning  techniques,  which  essentially 
attempt  to  convert  the  system  into  a  better  conditioned  matrix  equation  with  the 
same  solution  [71]  -  [73].  Still  another  technique,  the  reorthogonal izat ion  of 
the  p-functions  [70],  generally  requires  the  storage  of  the  set  {p  ^ }  and  thus  is 
not  immediately  compatible  with  the  goals  of  this  study. 

3.3.  TM  Scattering  from  Perfectly  Conducting  Cylinders 

One  of  the  simplest  moment  method  formulations  is  the  scattering  of  TM 
plane  waves  from  an  infinite,  perfectly  conducting  cylinder.  The  formulation  is 
general  and  can  treat  parallel  cylinders  of  arbitrary  shape,  including  thin 
strips.  The  matrix  equation  arising  from  a  moment  method  discretization  of  the 
electric-field  integral  equation  (EFIE)  has  been  discussed  by  Harrington  [74]. 

An  implementation  of  the  CGM  for  the  solution  of  the  matrix  equation  and  a 
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discussion  of  the  accuracy  of  the  numerical  solution  are  described  in  a  report 
by  Peterson  and  Mittra  [57],  and  briefly  summarized  in  Chapter  4. 

The  following  figures  illustrate  the  rate  of  convergence  for  several 
examples  of  the  above  moment  method  formulation.  Since  in  theory  the  CGM  and  OGM 
terminate  in  N  steps  or  less  for  an  order-N  system,  it  is  convenient  to  nor¬ 
malize  the  number  of  iterations  to  N.  This  approach  facilitates  the  observation 
of  trends  which  might  otherwise  be  obscured  by  examples  of  different  order. 

Thus,  the  figures  show  plots  of  the  residual  norm  as  defined  in  Equation 
(3.1)  versus  the  normalized  iteration  step  n/N. 

In  all  cases,  we  arbitrarily  consider  the  system  "solved"  once  a  solution 
is  produced  which  satisfies  <  0.0001.  Note  that  in  practice,  if  the  matrix 
involved  is  ill-conditioned,  the  residual  norm  is  not  always  a  valid  indication 
of  an  accurate  solution.  Under  these  conditions,  it  may  be  advisable  to  ter¬ 
minate  the  algorithm  at  a  lower  value  of  N  . 

6  n 

Various  parameters  describing  the  data  presented  in  these  figures  have  been 
controlled  in  order  to  investigate  any  general  effects  they  may  have  on  relative 
rates  of  convergence.  These  parameters  include  the  cell  density  per  wavelength 
used  within  the  moment  method  discretization,  the  presence  or  absence  of  simple 
symmetries  in  the  equation  and  solution,  and  the  degree  to  which  different  sized 
cells  (or  in  later  sections,  mixtures  of  different  permittivities  representing 
inhomogeneous  dielectric  materials)  are  incorporated  into  the  model. 

Furthermore,  we  consider  examples  of  the  CGM,  OGM,  and  AIGM.  The  specific  form 
of  the  AIGM  in  use  has  been  discussed  in  Section  3.1. 

Figure  3.1  shows  the  convergence  of  the  CGM  for  three  examples,  all  with 
the  cell  density  fixed  at  10.0  cells/XQ.  Two  of  these  systems  possessed  a  degree 
of  symmetry,  resulting  from  symmetry  planes  in  both  the  scatterer  geometry  and 
in  the  excitation  being  preserved  by  the  discretization.  Often,  the  CGM  is 
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observed  Co  converge  at  a  faster  rate  for  systems  with  this  type  of  symmetry  as 
compared  to  similar  equations  with  no  symmetry.  An  examination  of  the  expansion 
functions  produced  by  the  CGM  for  a  symmetric  example  showed  that  the  p- 
functions  contained  the  same  type  of  symmetry  as  the  solution.  This  may  explain 
the  quicker  convergence,  since  only  half  of  the  normal  set  of  expansion  func¬ 
tions  is  needed  to  span  the  symmetric  part  of  the  solution  space. 

Figure  3.2  shows  the  convergence  of  the  CGM  for  three  examples  of  moment 
method  systems  created  with  different  cell  densities.  These  equations  all 
posessed  a  degree  of  symmetry.  The  data  illustrate  a  general  trend  which  has 
been  observed  in  a  wide  variety  of  examples,  namely,  that  the  rate  of  convergence 
is  dependent  upon  the  sampling  density,  and  is  usually  faster  for  problems  with 
a  higher  density  of  cells  per  free-space  wavelength. 

The  interdependence  of  convergence  rate  and  sampling  density  is  related  to 
the  discussion  in  Section  3.2  concerning  the  eigenvalues  of  the  integral  opera¬ 
tor  and  their  projections  into  the  matrix  operator.  The  concept  is  better 
illustrated  by  Table  3.1,  which  shows  values  of  the  residual  norm  produced  by 
the  CGM  for  a  system  representing  a  circular  cylinder  illuminated  by  a  plane 
wave.  Four  different  levels  of  discretization  are  shown.  Note  that  while  only 
three  iterations  are  required  to  solve  the  4x4  system,  5  iterations  are  required 
to  solve  the  8x8,  16x16,  and  32x32  systems.  This  process  was  carried  out  for 
84x64  and  128x128  matrix  representations  as  well  (not  shown)  for  which  the  CGM 
also  required  5  iterations  to  reduce  the  residual  norm  below  0.0001.  The  impli¬ 
cations  here  are  that  the  corresponding  integral  operator  has  five  independent 
eigenvalues  and  that  the  matrix  is  "converging"  toward  the  integral  as  the 
order  of  the  system  increases.  (This  particular  geometry  yields  an  exact  solu¬ 
tion,  and  the  numerical  solution  appears  to  converge  to  the  exact  result  as  the 
order  of  the  system  increased.  The  comparison  is  presented  in  Section  4.3.)  Note 


41 


TABLE  3.1 

VALUES  OF  THE  RESIDUAL  NORM  PRODUCED  BY  THE  CGM  VERSUS  ITERATION  STEP 
FOR  FOUR  DIFFERENT  DISCRETIZATIONS  OF  THE  SAME  INTEGRAL  EQUATION 


order  of  system 

N  =  4 

8 

16 

32 

cell  density 

4.0  “US 

A0 

8.0 

16.0 

32.0 

a  =  0 

N  -1.0 

n 

1.0 

1.0 

1.0 

i 

0.359 

0.366 

0.361 

0.358 

2 

0.100 

0.114 

0.115 

0.115 

3 

8.9  x  10'10 

0.0142 

0.0161 

0.0161 

4 

- 

0.000706 

0.00128 

0.00132 

5 

- 

2.2  x  10~7 

_ 

6.9  x  l(f5 

8.0  x  10-5 

42 


also  that  the  value  of  the  residual  norm  at  each  step  appears  to  converge  as  the 
cell  density  increases,  possibly  to  the  value  of  the  residual  norm  which  would 
be  produced  by  an  application  of  the  CGM  to  the  integral  equation  itself  (if 
such  an  application  could  be  carried  out  analytically,  which  has  not  been  done 
to  date).  In  practice,  a  cell  density  of  10.0  cells/XQ  is  usually  considered 
necessary  to  produce  reasonably  accurate  numerical  solutions  to  these  integral 
equations.  According  to  Table  3.1,  the  numerical  values  of  the  residual  norm 
stabilize  at  about  this  cell  density.  Additional  examples  of  the  type  shown  in 
Table  3.1  are  presented  in  later  sections. 

Figure  3.3  shows  the  convergence  rates  of  2  examples  constructed  with  mixed 
cell  densities  in  the  "accurate"  range.  The  particular  models  in  use  involved 
large  and  small  cells  placed  immediately  adjacent  to  each  other,  and  contained 
no  simple  symmetry  planes.  Note  that  the  rate  of  convergence  in  these  examples 
is  slower  than  in  the  previous  examples.  Jennings  [30]  suggests  that  the  con¬ 
vergence  of  the  CGM  can  be  improved  by  scaling  so  that  the  main  diagonal  of  the 
matrix  has  identical  elements.  No  attempt  was  made  here  to  incorporate  scaling, 
which  may  explain  the  relatively  slower  convergence  for  the  mixed  cell  cases. 

Figure  3.4  shows  the  convergence  behavior  for  2  other  symmetric  systems,  one 
corresponding  to  a  circular  cylinder  and  the  other  a  pie-shaped  cylinder.  Both 
of  these  scatterers  are  also  internal  Ly  resonant  cavities,  such  that  the 
integral  equation  describing  the  external  scattering  problem  fails  because  of 
the  nontrivial  homogeneous  solution  due  to  the  internal  cavity  problem.  In  other 
words,  there  is  no  unique  solution  to  the  EF1E  for  the  specific  geometries  under 
consideration.  Although  the  CGM  was  able  to  accurately  solve  the  matrix 
equations  in  these  cases,  (as  was  verified  by  independent  computations  and  not 
merely  the  value  of  the  residual  norm),  the  numerical  solution  was  not  the 
desired  result  for  the  external  scattering  problem.  The  "internal  resonance" 


Figure  3.3. 


Convergence  of  the  CGM  for  examples  of  PEC  cylinders  TM 
polarization.  ~ ' 

N  *  63,  10.0-50.0  cells/Xg,  no  symmetry,  1  bent  strip, 

1  flat  strip,  mixed  cell  sizes  on  each  strip 
N  =*  51,  10.0-50.0  cells/Xg,  no  symmetry,  2  coplanar  strips 


Convergence  of  Che  CGM  for  examples  of  PEC  cylinders,  TM 
polarization,  where  geometries  are  "internally  resonant.' 
—  ■  •  -  N  *  67,  3.4  cells/Ag,  symmetry,  circular  cylinder 
-  N  =  120,  10. 0  cells/Xg*  symmetry,  pie-shaped  cylinde 
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situation  and  various  remedies  have  been  discussed  elsewhere  [75)  -  [78].  In 
practice,  the  characteristic  leveling  of  the  curves  shown  in  Figure  3.4  indicates 
an  ill-conditioning  problem,  and  this  behavior  may  be  useful  in  some  instances 
as  a  flag  to  identify  problem  situations.  In  spite  of  the  ill-conditioning, 
round-off  errors  do  not  appear  problematic  in  these  examples,  which  did  converge 
at  relatively  quick  rates.  The  symmetry  in  these  systems  may  be  responsible  for 
the  quick  convergence. 

Figures  3.5  and  3.6  show  the  performance  of  the  OGM  for  seve'al  examples. 

As  expected,  the  residual  norm  does  not  always  decrease  monotonical Iv  with  the 
OGM,  and  rates  of  convergence  are  similar  to  the  CGM.  In  fact,  the  CGM  and  OGM 
usually  require  the  same  number  of  iterations  to  solve  a  given  system,  which  is 
to  be  expected  because  of  the  "finite-step  convergence  in  a  number  of  steps 
equal  to  the  number  of  eigenvalues"  property  shared  by  both  of  these  algorithms. 
Figure  3.7  shows  the  residual  norm  produced  by  the  CGM  and  OGM  for  a  given 
example  involving  highly  mixed  cell  sizes  and  no  symmetry  planes.  Both 
algorithms  converge  at  similar  rates. 

Figure  3.8  shows  convergence  data  obtained  with  the  AIGM.  As  implemented 
here,  the  AIGM  is  observed  to  converge  at  a  slower  rate  than  the  CGM  and  OGM 
algorithms.  Figure  3.9  shows  a  comparison  of  the  CGM  and  AIGM  for  the  same 
system.  Figure  3.10  compares  all  three  algorithms  for  an  ill-conditioned  system 
corresponding  to  an  "internally  resonant"  circular  cylinder.  The  convergence  of 
the  AIGM  in  this  case  is  very  slow.  The  behavior  of  the  OGM  illustrated  in 
Figure  3.10  is  interesting  because  the  residual  norm  decreases,  then  increases 
to  a  value  exceeding  its  initial  value.  Although  this  behavior  is  observed  for 
a  special  type  of  example,  in  particular  one  where  the  condition  number  of  the 
matrix  lAl  is  in  excess  of  300000,  it  suggests  that  it  may  be  misleading  to  ter¬ 
minate  the  OGM  based  upon  the  value  of  an  indicator  which  does  not  decrease 


Figure  3.5.  Convergence  of  the  OGM  for  examples  of  PEC  cylinders,  TM 
polarizat ion . 

"  N  ■  85,  10.4  cells/A.^,  symmetry,  circular  cylinder 
N  *  30,  9.5  cells/Xg,  symmetry,  circular  cylinder 

-  N  *  64,  10.0  cells/Ag,  no  symmetry,  semi-circular  cylinder 
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Convergence  of  the  OGM  for  examples  of  PEC  cylinders,  TM 
po  lar izat ion . 

-  N  *  75,  10.0  cells/Xg,  symmetry,  circular  cylinder 

-  N  *  64,  10.0  cells/Xg,  no  symmetry,  semi-circular 

cylinder  with  2  coplanar  strips 
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Figure  3.8.  Convergence  of  the  AIGM  for  examples  of  PEC  cylinders,  TM 
polar izat  ion . 

N  *  60,  10.0  cells/Ag,  no  symmetry,  2  coplanar  strips 

-  N  3  60,  100.0  cells/Ag,  symmetry,  2  coplanar  strips 

-  N  3  20,  10. 0  cells/Ag,  no  symmetry,  2  coplanar  strips 


_n 

N 


Figure  3.10.  Comparison  of  the  CGM,  OGM ,  and  AIGM  for  an  example  of  a  PEC 

cylinder,  TM  polarization,  where  geometry  is  "internally  resonant 
circular  cylinder,  N  ■  40,  7.8  cells/Ag,  symmetry. 
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monotonically.  It  also  raises  some  questions  about  the  validity  of  the  residual 
norm  itself  as  an  indicator. 

3.4.  TM  Scattering  from  Dielectric  Cylinders 

A  scheme  for  the  analysis  of  arbitrary,  lossy  or  lossless  dielectric  cylin¬ 
ders  illuminated  with  plane  TM  waves  was  developed  by  Richmond  [79].  Here,  we 
use  the  Richmond  formulation  but  solve  the  resulting  moment  method  matrix 
equation  with  the  iterative  algorithms  of  Chapter  2. 

Figure  3.11  shows  the  rates  of  convergence  of  the  CGM  for  three  systems 
representing  dielectric  cylinders  discretized  with  different  cell  densities.  As 
was  observed  for  previous  examples,  a  general  trend  indicating  faster  con¬ 
vergence  when  smaller  cell  sizes  are  used  is  indicated.  For  cylinder  examples 
(infinite  cylinders  described  by  their  cross-sections)  involving  dielectric 
material,  cell  densities  refer  to  the  number  of  cells  per  unit  cross-sectional 
area,  where  area  vs  measured  in  terras  of  the  wavelength  in  the  dielectric 
material.  If  \Q  is  the  free-space  wavelength,  is  the  wavelength  in  the 
dielectric  material  defined  by 

xo  =  ^TTT  xd  o.io) 

where  is  the  complex  relative  permittivity.  Note  that  the  cell  density  of 
2 

9.6  cells/A^  used  in  one  of  the  examples  is  normally  considered  far  below  what 
is  required  for  the  moment-method  discretization  to  accurately  model  the  scat¬ 
tering  problem. 

Figure  3.12  shows  data  for  2  examples  where  mixed  permittivities  were  used 


in  the  model.  This  did  not  appear  to  slow  the  rate  of  convergence  over  cases 
where  a  constant  permittivity  was  used  throughout  the  model. 
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Figure  3.11.  Convergence  of  Che  CGM  for  examples  of  dielectric  cylinders,  TM 
po l ar izat ion . 

-  N  =  63,  835  cells/A^,  homogeneous,  rectangular  cylinder, 

ef  =  3 

N  =  21,  104  cells/A^,  homogeneous,  circular  cylinder, 

£r  =  2.56 

-  N  =  25,  9.6  cells/A^,  homogeneous,  square  cylinder, 

=  2.56 
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Figures  3.13  and  3.14  show  the  performance  of  the  OGM,  CGM,  and  AIGM 
applied  to  geometries  modeled  with  mixed  cell  sizes  and  permittivities.  Again, 
previous  remarks  concerning  the  algorithms  apply. 

Table  3.2  shows  the  residual  norm  versus  the  unnormalized  iteration  step 
for  three  different  models  of  a  homogeneous,  circular  dielectric  cylinder,  illu¬ 
minated  by  a  plane  TM  wave.  The  three  models  consist  of  21,  61,  and  101  square 
celts  configured  to  approximate  the  circular  cross-section  of  the  cylinder.  As 
the  number  of  cells  in  the  model  increases,  the  order  of  the  system  increases  as 
well,  but  in  each  case  only  5  iterations  are  required  for  CGM  solution.  As 
discussed  in  Section  3.3,  the  matrix  representation  appears  to  be  "converging" 
to  something  as  the  discretization  is  refined.  This  example  is  studied  in 
Chapter  4,  where  it  is  shown  that  the  solution  to  the  discrete  system  appears  to 
converge  to  the  analytical  solution  for  the  circular  cylinder. 

3 . 3  TE-wave  Scattering  from  Dielectric  Cylinders 

A  computational  method  for  TE-wave  scattering  from  arbitrary  dielectric 
cylinders  was  formulated  by  Richmond  [80].  This  approach  uses  the  EFIE  with 
pulse  basis  functions  and  point  matching,  and  an  iterative  implementation  of  this 
scheme  is  developed  in  Chapter  6.  In  this  section,  the  CGM  and  AIGM  are  applied 
to  solve  the  moment  method  system,  and  their  rates  of  convergence  are  presented 
via  graphs  showing  the  residual  norm  versus  the  normalized  iteration  step.  In 
general,  these  rates  of  convergence  depend  upon  the  cell  density  of  the  par¬ 
ticular  model  as  did  the  previous  examples  of  this  chapter.  The  cell  density  is 
defined  in  Section  3.4.  Note  that  the  convergence  rates  of  the  AIGM  are  slightly 
faster  than  the  CGM  for  these  equations.  As  throughout  the  rest  of  the  chapter, 
the  AIGM  is  implemented  using  an  approximate  inverse  obtained  by  inverting  the 
main  diagonal  of  the  matrix. 
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Figure  3  13.  Convergence  of  Che  CGM,  OGM,  and  AIGM  for  examples  of  dielectric 
cylinders,  TM  polarization. 

N  *  78,  225-645  cells/X^,  inhomogeneous  cylinder  with  er  ranging 
from  2—  j 1  to  4.5-j2 

-  CGM  and  OGM 

— —  AIGM 
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Figure  3.15  shows  the  convergence  rates  of  the  CGM  when  applied  to  several 
examples,  each  modeled  with  similar  cell  densities.  These  examples  serve  to 
emphasize  the  fact  that  while  we  observe  general  trends,  individual  examples  can 
deviate  somewhat  from  these. 

Figure  3.16  shows  convergence  rates  of  the  CGM  for  three  examples,  each 

with  different  cell  densities.  These  examples  again  illustrate  the  general  trend 

that  convergence  is  usually  faster  for  models  having  higher  cell  densities.  Note 

2 

that  the  cell  density  of  18.5  cells/A^  is  considered  inadequate  to  produce 
accurate  numerical  solutions  to  the  scattering  problem. 

Table  3.3  shows  values  of  the  residual  norm  versus  the  unnormalized  itera¬ 
tion  step  for  three  models  of  the  same  scattering  problem,  in  this  case  a  homo¬ 
geneous  circular  cylinder.  As  the  cell  density  increases,  the  number  of 
iterations  required  stabilizes,  and  the  values  of  the  residual  norm  n  appear  to 
stabilize  also.  As  discussed  in  Section  3.3,  we  interpret  this  as  an  indication 
that  the  discrete  system  is  "converging"  as  the  order  of  the  d isc ret izat ion  is 
refined.  The  question  of  whether  the  discrete  system  is  converging  to  the 
desired  integral  equation  is  reserved  for  discussion  in  Chapter  6.  Observe  that 

the  process  seems  to  stabilize  once  the  cell  density  exceeds  about 
2 

100  cells/A^,  an  indication  that  the  eigenvalue  structure  of  the  matrix  has 
stabil ized . 

Figure  3.17  shows  the  convergence  of  the  CGM  when  applied  to  models  con¬ 
sisting  of  mixed  cell  sizes  and  permittivities,  and  posessing  no  simple  types  of 
symmetry . 

Figure  3.18  shows  the  convergence  rates  of  the  A1GM  when  applied  to  several 
models  of  dielectric  cylinders.  Figures  3.19  and  3.20  compare  the  convergence 
rates  of  the  CGM  and  AIGM  for  two  matrix  equations.  These  figures  illustrate  a 
trend  which  has  been  observed  for  many  examples  of  the  TE  dielectric  cylinder 


Figure  3.15.  Convergence  of  the  CGM  for  examples  of  dielectric  cylinders, 

IE  polarization. 

N  *  168,  104  cells/X$,  homogeneous  circular  cylinder  with 
er  *  2.56 

-  N  ■  80,  132  cells/X^,  homogeneous  rectangular  cylinder  with 

er  *  3-j0.2 

-  N  *  42,  104  cells/Xj,  homogeneous  circular  cylinder  with 


TABLE  3.3. 


CONVERGENCE  OF  THE  CGM  FOR  THREE  MODELS  OF  A  CIRCULAR  DIELECTRIC 
CYLINDER  WITH  RADIUS  =  0.3183  A  ,  e  =  2.56,  TE  POLARIZATION 


order  of  system 

N  =  42 

168 

' 

672 

cell  density 

cells 

‘2 

104 

412 

n  31  0 

N  =  1 

1 

> 

i 

0.587 

0.585 

0.584 

2 

0.320 

0.340 

0.347 

3 

0.249 

0.258 

0.263 

4 

0.156 

0.178 

0.181 

5 

0.0881 

0.1 18 

0.124 

6 

0.0564 

0.0776 

0.0831 

7 

0.0372 

0.0513 

0.0570 

8 

0.0250 

0.0366 

0.0412 

9 

0.0127 

0.0235 

0.0272 

10 

0.00723 

0.0146 

0.0166 

11 

0.00359 

0.00880 

0.0107 

12 

0.00214 

0.00497 

0.00652 

13 

0.00137 

0.00270 

0.00361 

14 

0.000649 

0.00161 

0.00191 

15 

0.000238 

0.000941 

0.00108 

16 

0.0000876 

0.000638 

0.000663 

17 

- 

0.000378 

0.000457 

18 

- 

0.000171 

0.000258 

Convergence  of  the  CGM  for  examples  of  dielectric  cylinders, 

TE  polarization. 

N  *  90,  125-1063  cells/X^,  5  square  dielectric  cylinders 
in  close  proximity,  £r  *  2  -  jO  to  er  3  5  -  j0.8. 

-  N  3  80,  100-196  cells/X^,  inhomogeneous  skewed-rectangular 

cylinder  with  €r  3  2-j0.4,  3-j0.2,  and  4-jO.l 

-  N  3  80,  992-1962  cells/X^,  inhoraog  eneous  skewed-rectangular 

cylinder  with  er  3  5-jl  and  10-jl 


Convergence  of  the  AIGM  for  examples  of  dielectric  cylinders,  TE 
polar  izat ion . 

N  *  90,  365  cells/X^,  homogeneous  rectangular  cylinder  with 
er  *  3- jO . 5 

-  N  =  42,  100  cells/A^,  homogeneous  circular  cylinder  with 

er  3  5-j2 

-  N  =  42,  98  cells/Ag,  homogeneous  circular  cylinder  with 

er  3  5- j 5 
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problem,  namely,  that  the  A1GM  appears  to  converge  at  a  faster  rate  than  the  CGM. 
The  convergence  of  the  CGM  as  illustrated  in  Figure  3.20  is  much  slower  than 
would  normally  be  expected,  considering  the  previous  performance  of  the  CGM  for 
similar  problems.  In  fact,  the  plot  of  Figure  3.20  bears  some  resemblance  to  the 
plots  of  the  "internal  resonance"  problems  of  Figures  3.4  and  3.10,  suggesting 
that  there  might  be  a  problem  with  the  accuracy  of  the  numerical  result  in  this 
case . 

To  facilitate  a  study  of  this  type  of  behavior,  Figure  3.21  shows  the  per¬ 
formances  of  the  CGM  and  AIGM  for  an  example  where  the  cell  density  is  clearly 
inadequate  to  represent  the  scattering  problem  in  question.  As  should  be 
expected,  the  convergence  of  both  algorithms  is  very  slow.  The  slow  convergence 
of  the  CGM  is  indicative  of  inadequate  sampling  of  the  eigenvalue  structure  of 
the  integral  equation;  yet  the  process  converges  in  this  case  because  round-off 
errors  are  apparently  not  severe  enough  to  prevent  the  finite-step  termination 
of  the  algorithm.  The  behavior  of  the  CGM  is  very  similar  in  the  examples  of 
Figures  3.20  and  3.21,  yet  the  AIGM  behavior  is  very  different.  This  is  probably 
related  to  the  accuracy  of  the  approximate  inverse  in  the  case  where  |  er  |  becomes 
large.  In  both  cases,  the  performance  of  the  CGM  suggests  that  there  might 
be  a  problem  with  the  matrix  equation  in  use.  Further  study  has  shown  that  the 
accuracy  of  the  partLCular  moment  method  formulation  is  poor  for  large  values  of 
|--rl,  and  this  topic  is  examined  in  detail  in  Chapter  6.  Although  we  have  no 
reason  a  priori  to  expect  the  formulation  to  fail  in  this  case,  as  it  turns  out 
there  is  a  problem  which  might  have  gone  unnoticed  were  we  not  familiar  with  the 
normal  convergence  rate  of  the  CGM. 

The  convergence  rate  of  the  CGM  for  an  additional  example  of  a  problem 
which  was  expected  to  adequately  represent  a  scattering  problem  is  depicted 
in  Figure  3.22.  The  numerical  solution  to  the  matrix  in  question  bore  no 
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resemblence  to  the  analytical  solution,  which  was  available  for  this  geometry  (a 
homogeneous,  circular  cylinder  illuminated  by  a  plane  TE  wave).  The  permittivity 
involved  here  was  =  76-j278,  in  the  range  often  used  to  model  certain 
biological  media.  More  information  on  the  behavior  of  this  problem  is  presented 
in  Chapter  6. 

3.6.  Scattering  from  Cylinders  Containing  Both  Perfectly 

Conducting  and  Dielectric  Materials  for  Both  TM  and  TE 

Polarizations 

Chapter  4  considers  the  iterative  implementation  of  a  moment  method  analy¬ 
sis  for  cylinders  modeled  by  both  perfectly  conducting  (PEC)  and  dielectric 
materials.  Specifically,  the  TM  polarization  is  treated  using  the  EF1E  with 
pulse  basis  functions  and  point  matching.  The  TE  polarization  is  treated  by  com¬ 
bining  the  EFIE  for  the  dielectric  material  with  the  magnetic  field  integral 
equation  (MF1E)  for  the  PEC  material,  and  again  using  pulse  basis  functions  and 
point  matching.  Because  the  resulting  matrix  equations  contain  terms  repre¬ 
senting  the  interaction  between  the  different  types  of  material  present,  they 
are  somewhat  more  complicated  than  the  systems  presented  thus  far  in  this 
chapter.  For  instance,  the  matrix  equations  considered  in  Sections  3,3  and  3.4 
were  of  the  form 
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while  the  more  complicated  example  of  Section  3.5  involved  a  system  of  the  form 
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Equation  (3.11)  represents  a  scalar  integral  equation.  Equation  (3.12)  repre¬ 
sents  a  vector  equation  which  has  been  separated  into  x  and  y  components.  The 

blocks  labelled  G  and  G  represent  interactions  between  different  vector  com- 
xy  yx 

ponents.  In  the  combined  dielectric  and  PEC  case  the  matrix  equation  takes  the 
form 


(3.13) 
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for  the  TE  polarization.  The  sub-matrices  appearing  in  Equations  (3.13)  and 
(3.14)  represent  interactions  between  different  vector  components  or  different 
materials,  and  are  generally  different  orders  of  magnitude  than  the  blocks 
located  aLong  the  main  diagonal.  This  affects  the  condition  of  the  equation, 
and  prompts  us  to  consider  scaling  the  different  parts  of  the  matrix  to  improve 
the  convergence  of  the  iterative  algorithms  when  used  to  solve  these  systems. 
This  section  shows  the  performance  of  the  CGM  and  AIGM  when  applied  to  Equations 
(3.13)  and  (3.14).  Additional  information  on  the  discretizations  in  use  may  be 


found  in  Chapter  4. 

Consider  the  TM  polarization  described  by  Equation  (3.13).  In  order  to 
alleviate  any  problems  which  may  arise  due  to  the  presence  of  different  orders- 
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of-magnitude  in  the  blocks  which  make  up  the  matrix,  a  scale  factor  is  intro¬ 
duced  to  normalize  the  different  sub-matrices  to  the  same  order  of  magnitude. 
Equation  (3.13)  takes  the  form 


(3.15) 


z 


(3.16) 


and  the  scale  factor  S  is  chosen  to  give  the  two  types  of  unknowns  and  J  , 

.  TM 

the  same  order  of  magnitude.  For  the  following  examples,  the  value  S  -  0.025 

was  used  with  the  specific  form  of  the  equation  described  in  Chapter  4. 

Figure  3.23  shows  the  convergence  of  the  CGM  for  three  different  circular 
PEC  cylinders  coated  in  each  case  with  a  single,  uniform  layer  of  lossless 
dielectric  material.  Overall,  the  convergence  behavior  is  similar  to  that  seen 
in  previous  sections  of  this  chapter.  Note  that  in  keeping  with  the  previous 
results  there  is  a  general  trend  of  relatively  quick  convergence  (normalized  to 
the  order  of  the  system)  as  higher  cell  densities  are  used  in  the  models. 

Figure  3.24  shows  the  convergence  of  the  CGM  for  three  examples  involving  a 
single,  circular  PEC  cylinder  in  the  presence  of  one  or  more  rectangular 
dielectric  cylinders. 

Figure  3.25  compares  the  convergence  of  the  CGM  and  AIGM  for  a  matrix 
equation  representing  a  circular  PEC  cylinder  with  a  uniform  dielectric  coaling. 

Although  we  introduced  a  scale  factor,  experimentation  showed  that  the  rate 
of  convergence  of  the  iterative  algorithms  was  not  a  strong  function  of  the  spe¬ 
cific  scale  factor  in  use,  for  the  TM  polarization.  The  system  representing  the 


_n 

N 

Figure  3.23.  Convergence  of  the  CGM  for  examples  of  combination  dielectric/PEC 

cylinders,  TM  polarization,  circular  PEC  cylinders  with  homogeneous 
dielectric  coating,  scale  factor  ■*  0.025. 

-  N  *  99,  66.0  cells/Xg  PEC,  6180  cells/X^  dielectric 

N  -  99,  16.5  cells/Xg  PEC,  384  cells/Xjj  dielectric 
-  N  *  38,  10.0  cells/ Xg  PEC,  151  cells/Xs  dielectric 


Convergence  of  the  CGM  for  examples  of  combination  dielectric/PEC 
cylinders,  TM  polarization,  with  scale  factor  *  0.025. 

N  *  272,  10.0  cells/XQ  PEC,  110  cells/X^  dielectric,  circular 
PEC  cylinder  behind  square  cylinder,  with  er  *  2  —  j 0 . 3 

-  N  *  60,  10.0  cells/XQ  PEC,  132-279  cells/X^  dielectric, 

circular  PEC  cylinder  near  3  square  dielectric  cylinders  with 
tr  -  1.5-j0. 1,  3.0-j0. 3,  4-j0 

— .  N  *  52,  10.0  cells/XQ  PEC,  269-424  cells/X^  dielectric, 

circular  PEC  cylinder  near  2  rectangular  dielectric  cylinders 
with  er  *  3-J0.5,  4-j0.8 
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TE  polarization  appears  to  be  much  more  sensitive  to  proper  scaling,  as 
illustrated  below. 

Equation  (3.14),  representing  the  TE  polarization,  can  be  scaled  to  take 
the  form 


Jp  =  s’p'  p  (3.18) 

TE  TE 

The  scale  factors  =  0.01  and  S2  =  377  appear  to  work  well  with  the  CGM  to 
ensure  normal  convergence  behavior.  As  an  example,  Figure  3.26  shows  the  con¬ 
vergence  of  the  CGM  for  a  circular  PEC  cylinder  with  a  uniform  dielectric 
coating.  Note  that  without  scaling,  the  CGM  does  not  converge  to  a  solution. 

Figure  3.27  shows  the  convergence  of  the  CGM  for  geometries  involving  a 
single  circular  PEC  cylinder  in  the  presence  of  one  or  more  rectangular 
dielectric  cylinders.  Typically,  convergence  was  fairly  rapid  for  this  type  of 
mode  1  . 

Figure  3.28  shows  the  convergence  of  the  CGM  and  AIGM  for  an  example 
involving  a  circular  PEC  cylinder  and  two  square  dielectric  cylinders.  Although 
the  AIGM  is  unable  to  produce  a  solution  with  V  <  0.0001  in  this  case,  the 
solution  found  after  36  iterations  was  reasonably  accurate.  The  slow  convergence 
of  the  CGM  for  this  example  suggests  the  possibility  of  a  problem;  in  fact,  the 
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N 


Figure  3.26.  Convergence  of  CGM  for  an  example  of  a  circular  PEC  cylinder  coated 
with  a  homogeneous  dielectric  layer,  TE  polarization. 

N  -  48,  25.5  cel  Is/ Xq  PEC,  137  cells/X§  dielectric,  er  =  3-jO 
-  Sp  -  0.01  SP  -  377 

- sp  -  1  sp  -  1 


Figure  3.28.  Comparison  of  the  convergence  of  the  CGM  and  AIGM  for  an  example  of 
a  cylinder  containing  both  dielectric  and  PEC  material.  N  =  90, 
10.0  celis/Xg  PEC,  204-369  cells/X^  dielectric,  £r  a  3-j0  and 
3- jO .  5 ,  Sp  =■  0.01,  Sp  “  377  ,  circular  PEC  cylinder  and  2 
square  dielectric  cylinders. 
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performance  of  this  numerical  approach  is  very  sensitive  to  some  of  the  approxi 
mations  made  in  the  evaluation  of  the  matrix  elements.  A  detailed  discussion  of 
these  approximat ions  and  their  effects  on  the  accuracy  of  the  final  solution 
is  provided  in  Chapter  4.  The  particular  example  illustrated  in  Figure  3.28  is 
borderline  in  that  both  the  rate  of  convergence  and  the  accuracy  of  the  numeri¬ 
cal  solution  could  be  improved  somewhat  by  better  evaluation  of  the  matrix 
e lements . 

3.7.  Use  of  the  CGM  to  Treat  Multiple  Excitations 

As  mentioned  in  Section  3.1,  the  amount  of  calculation  per  iteration  per¬ 
formed  by  the  CGM  algorithm  is  roughly  divided  between  the  tasks  of  generating 
the  expansion  functions  (the  p-functions  of  Section  2.3)  and  computing  their 
coefficients.  Thus,  one  way  of  efficiently  treating  multiple  right-hand  sides 
is  to  simultaneously  expand  several  solutions  in  terms  of  a  single  set  of  expan 
sion  functions.  In  theory,  this  approach  will  permit  us  to  save  about  half  of 
the  cost  per  iteration  of  treating  each  additional  solution  ( corresponu ing  to 
each  additional  right-hand  side)  without  adding  significantly  to  the  storage 
requirements.  In  order  to  test  this  idea  in  practice,  the  following  data  were 
generated . 

Figure  3.29  shows  the  rate  of  convergence  of  the  CGM  as  illustrated  by  the 
residual  norm  for  five  different  systems  arising  in  the  analysis  of  a  single, 
circular  PEC  cylinder  by  the  moment  method  approach  of  Section  3.3.  Each  solu¬ 
tion  corresponds  to  a  different  excitation,  in  this  case,  plane  waves  impinging 
on  the  scatterer  from  five  different  angles.  The  p-functions  from  the  0  degree 
incidence  wave  were  used  to  expand  all  five  solutions.  The  process,  however, 
evidently  does  not  vield  the  desired  improvement  in  efficiency,  as  only  the  0 
degree  solution  is  obtained  to  the  necessary  accuracy.  Some  of  the  other 
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Figure  3.29.  Convergence  of  the  CGM  modified  to  treat  simultaneous  excitations 
for  an  example  of  a  circular  PEC  cylinder,  TM  polarizat ion .  N  * 

9.6  cells/XQ* 

-  9i  ,  o* 

.  91  -  2.5° 
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solutions  are  clearly  no  better  than  their  initial  estimate,  even  after  the 
solution  to  0  degree  incidence  is  obtained. 

Figure  3.30,  representing  to  a  different  scatterer,  illustrates  similar 
behavior  when  the  CGM  is  applied  to  the  simultaneous  treatment  of  three  dif¬ 
ferent  excitations.  Clearly,  in  this  case  it  would  be  more  efficient  to  treat 
each  system  separately. 

The  reasons  for  the  failure  of  the  above  technique  are  necessarily  related 
to  the  strong  point  of  the  CGM  —  namely,  its  quick  convergence  in  most  cases. 
Because  the  expansion  functions  are  geared  to  represent  the  solutions 
corresponding  to  the  specific  right-hand  side  in  use,  the  convergence  of  the 
algorithm  is  relatively  fast.  These  functions,  however,  are  not  geared  to  repre¬ 
sent  any  function,  even  though  in  theory  a  finite  number  of  them  span  the  space. 
Because  of  round-off  errors  affecting  the  orthogonality,  they  apparently  do  not 
span  the  space  (in  spite  of  the  fact  that  they  can  represent  the  single  solution 
they  were  generated  for).  In  addition,  it  has  been  mentioned  that  the  p- 
functions  may  Lake  on  the  symmetries,  if  any,  of  the  corresponding  solutions. 
Thus,  they  clearly  would  not  be  useful  for  the  representation  of  solutions  which 
did  not  possess  the  same  symmetries.  For  all  of  these  reasons,  the  simultaneous 
expansion  of  several  solutions  does  not  appear  effective  for  the  efficient 
treatment  of  multiple  systems. 

An  alternate  approach,  that  of  using  a  previous  solution  to  generate  an 
initial  estimate  of  a  different  solution,  may  prove  feasible  for  the  treatment 
of  multiple  right-hand  sides.  However,  as  evidenced  by  the  theoretical  con¬ 
vergence  rates  given  in  Section  3.2  and  by  observations  made  in  practice,  the 
rate  of  convergence  is  relatively  independent  of  the  initial  guess.  Thus,  a 
significant  improvement  in  efficiency  is  only  possible  by  finding  a  fairly 


accurate  initial  estimate  of  the  solution. 
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Figure  3.30.  Convergence  of  the  CGM  modified  to  treat  simultaneous  excitations 
for  an  example  of  a  circular  PEC  cylinder,  TM  polarization. 

N  -  35,  9.3  cel  Is/ Xq . 
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Generally  speaking,  the  goals  of  fast  convergence  and  ability  to  treat 
multiple  right-hand  sides  are  somewhat  diametrically  opposed,  as  complete  flexi¬ 
bility  in  the  latter  mandates  that  the  full  set  of  expansion  functions  be 
available  (which  negates  the  former).  A  better  approach  for  dealing  with 
multiple  right-hand  sides  would  require  an  algorithm  which  does  not  converge 
fast  for  a  single  problem,  but  does  generate  useful  expansion  functions  for  the 
treatment  of  other  systems.  A  non-iterative  method  such  as  Gaussian  elimination 
may  be  the  best  currently  available  example  of  such  an  algorithm,  in  spite  of 
the  additional  storage  required  by  all  general  direct  solution  procedures. 

3.8.  Summary 

This  chapter  has  presented  a  discussion  of  the  theoretical  convergence  of 
the  CGM,  and  has  illustrated  the  convergence  of  all  three  of  the  iterative 
algorithms  of  Chapter  2  for  a  variety  of  electromagnetic  scattering  problems. 

The  examples  presented  were  selected  from  a  wide  range  of  test  cases,  and  are 
believed  to  represent  the  extremes  that  typically  arise  in  practice. 

Based  upon  these  data,  the  CGM  solution  of  a  moment  method  system  usually 
requires  N/4  to  N/2  iterations,  assuming  that  the  cell  densities  in  use  are  suf¬ 
ficient  to  ensure  adequate  sampling  of  the  original  integral  equation,  and  where 
N  is  the  order  of  the  matrix  equation.  However,  it  is  not  uncommon  for  the 
algorithm  to  require  as  many  as  3N/4  iterations,  especially  if  the  moment  method 
formulation  involves  mixed  cell  sizes  as  was  observed  in  Section  3.3.  As  higher 
cell  densities  are  used,  faster  rates  of  convergence  (relative  to  N)  are 
obtained.  If  convergence  is  much  slower  than  this,  or  if  the  residual  norm 
remains  virtually  constant  for  many  iterations,  it  is  possible  that  the  matrix 
is  not  an  accurate  model  of  the  original  equation  and  should  be  modified.  On  the 
other  hand,  very  fast  convergence,  i.e.,  N/10  iterations  or  less,  has  never  been 
observed  for  the  examples  considered  here  unless  the  cell  densities  in  use  are 
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extremely  high.  The  convergence  behavior  of  the  CGM  is  linked  to  the  eigenvalue 
structure  of  the  matrix,  which  itself  is  an  approximation  to  the  eigenvalue 
spectrum  of  the  integral  operator.  Thus,  the  convergence  behavior  of  the  CGM  is 
feedback  which  appears  to  be  useful  for  indicating  the  ad  nacy  of  the  discreti¬ 
zation  used  to  form  the  matrix.  These  conclusions  are  based  upon  observations 
made  with  several  types  of  integral  equations  representing  electromagnetic  scat¬ 
tering  problems  and  may  not  be  applicable  to  other  types  of  equations. 

The  AIGM  algorithm  converged  slower  than  the  CGM  and  OGM  for  all  of  the 
systems  except  the  TE  dielectric  cylinder  matrix  equation.  However,  this  example 
shows  that  even  with  the  simple  idea  of  generating  an  approximate  inverse  opera¬ 
tor  by  inverting  the  main  diagonal  of  the  matrix,  the  techniaue  is  feasible  and 
may  be  more  efficient  than  the  CGM.  It  may  be  advantageous  to  expend  additional 
effort  to  generate  a  better  inverse,  such  as  inverting  a  banded  approximation 
to  the  matrix.  For  a  specific  application,  more  information  can  be  brought  to 
bear  on  the  problem  of  finding  a  good  approximate  inverse.  Our  purpose  here  was 
to  examine  general  algorithms  that  were  not  tied  to  any  particular  geometry  or 
symmetry,  and  thus  we  did  not  attempt  to  find  a  better  approximate  inverse  for 
any  of  the  examples  considered. 
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4.  THE  MATRIX  ELEMENT  REGENERATION  (MER)  APPROACH 

4.1.  Introduction 

Iterative  techniques  have  been  employed  in  the  solution  of  electromagnetic 
scattering  problems  primarily  because  they  permit  any  sparseness  or  redundancy 
of  the  associated  matrix  equation  to  be  fully  exploited.  It  has  long  been  appre 
ciated  that  the  computer  storage  required  for  the  direct  solution  of  a  large 
system  places  a  practical  limit  on  the  electrical  size  of  the  scatterer  to  be 
analyzed.  Thus,  specialized  approaches  have  been  developed  which  build  sparse¬ 
ness  or  redundancy  into  the  matrix  equation,  permitting  more  efficient  use  of 
fast-access  memory.  Examples  of  these  approaches  are  discussed  in  Chapters  5  - 
7. 

The  discretizations  used  with  the  above-mentioned  "special"  approaches  limit 
the  scope  of  these  methods  in  most  cases  to  problems  involving  surfaces  of 
constant  curvature  or  volumes  which  can  be  represented  by  evenly-spaced  sub- 
domains.  The  present  chapter  investigates  a  different  implementation  of  itera¬ 
tive  algorithms,  one  that  is  suitable  for  the  treatment  of  arbitrary  geometries 
In  general,  these  methods  are  not  as  efficient  as  the  specialized  techniques  of 
Chapters  5-7,  because  they  require  more  computation  per  iteration  step. 
However,  they  may  provide  an  effective  alternative  for  geometries  which  are  not 
easily  treated  by  the  specialized  methods. 

4.2.  The  Matrix  Element  Regeneration  (MER)  Approach 

The  approach  to  be  investigated  is  based  upon  the  simple  idea  of 
"recomputing"  each  matrix  element  as  it  is  needed  in  the  iterative  solution  of 
the  matrix  equation.  We  denote  this  the  matrix  element  regeneration  (MER) 
approach.  The  goal  of  this  investigation  is  to  determine  if  the  MER  can  be 
implemented  efficiently  for  different  types  of  electromagnetic  scattering 
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problems.  Earlier  work  on  this  topic  dealt  primarily  with  problems  represented 
by  a  simple  scalar  integral  equation,  and  proved  encouraging  [54]  -  [57].  Here, 
the  MER  is  applied  to  a  more  complicated  integral  equation  representing  a  scat- 
terer  modeled  by  a  combination  of  perfectly  conducting  and  lossy  dielectric 
materials . 

The  idea  behind  the  MER  approach  is  that  it  may  be  more  efficient  to  generate 
each  matrix  element  as  it  is  needed  than  to  repeatedly  transfer  the  necessary 
numerical  values  from  an  out-of-core  storage  device  such  as  a  disk  or  magnetic 
tape  unit.  In  order  to  test  this  idea,  the  relative  execution  times  of  both 
approaches  are  compared. 

Since  simple  discretization  schemes,  such  as  the  moment  method  using  pulse 
basis  functions  and  point  matching,  are  likely  to  produce  matrix  equations  with 
relatively  simple  matrix  elements,  they  appear  to  be  the  best  candidates  for 
efficient  MER  implementation.  The  methods  to  be  considered  here  are  all  based 
upon  these  simple  basis  and  testing  functions.  However,  the  use  of  simple 
discretizations  has  sometimes  led  to  inaccuracy  in  the  final  numerical  solution 
[81],  In  order  to  determine  the  overall  value  of  the  methods,  the  accuracy  of 
each  numerical  approach  should  be  investigated  thoroughly  for  a  variety  of 
problems.  While  we  do  not  attempt  an  exhaustive  study  here,  several  examples 
are  given  in  order  to  compare  numerical  results  with  exact  analytical  solutions. 

The  objective  of  the  MER  process  is  to  reduce  the  amount  of  storage  required 
for  an  NxN  system  to  some  small  multiple  of  N.  Information  describing  the 
geometry  and  materials  can  be  stored  directly,  as  can  the  matrix  elements 
comprising  the  main  diagonal  of  the  NxN  system.  Savings  are  obtained  by  reducing 
the  required  storage  for  the  off-diagonal  matrix  elements. 

4.3.  TM-wave  Scattering  by  Conducting  Cylinders 

As  an  example,  consider  the  problem  of  TM-wave  scattering  by  perfectly  con¬ 
ducting  infinite  cylinders.  A  detailed  description  of  the  formulation  and 
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discretization  of  the  electric-field  integral  equation  are  provided  by 
Harrington  [74] .  Specifically,  the  approach  requires  a  conducting  surface  to  be 
modeled  by  subsectional  strips,  over  which  the  unknown  current  density  is 
treated  as  a  constant  (i.e.,  pulse  basis  functions).  The  equation  is  enforced  by 
point  matching  at  the  center  of  each  strip.  The  off-diagonal  matrix  elements  can 
be  found  approximately  by  replacing  the  pulse  basis  functions  by  Dirac  delta 
functions  located  at  the  center  of  each  strip,  and  take  the  form 
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In  Equation  (4.1),  represents  the  distance  from  the  center  of  strip  n  to  the 

2 

center  of  strip  m,  and  (*)  is  the  zeroth-order  Hankel  function  of  the  second 
kind.  The  Hankel  function  can  be  computed  efficiently  using  linear  interpolation 
from  a  look-up  table,  a  task  requiring  only  a  few  arithmetic  operations.  Note 
that  although  the  size  of  the  look-up  table  varies  in  proportion  with  the  maxi¬ 
mum  linear  dimension  of  the  scatterer  geometry  under  consideration,  it  is  typi¬ 
cally  much  smaller  than  the  NxN  matrix  it  replaces. 

Execution  time  data  for  small-order  systems  are  provided  in  a  report  by 
Peterson  and  Mittra  [57],  and  are  reproduced  in  Figure  4.1.  The  original  data 
show  the  execution  time  required  to  iteratively  solve  the  system  using  the  CGM 
and  to  compute  the  bistatic  radar  cross  section,  for  each  of  two  cases.  In  the 
first  case  the  entire  NxN  matrix  is  stored  in  fast-access  memory;  in  the  second 
case  all  of  the  off-diagonal  matrix  elements  are  recomputed  whenever  needed  in 
accordance  with  the  MER  procedure.  The  original  data  have  been  augmented  with 
data  showing  the  corresponding  execution  time  when  each  of  the  off-diagonal 
matrix  elements  is  transferred  from  a  disk  storage  device  as  needed.  The  CDC 
CYBER  175  computer  used  for  the  study  is  a  main-frame  time-sharing  machine,  and 
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the  data  may  be  machine  dependent  to  some  extent.  In  addition,  no  attempt  was 
made  to  transfer  blocks  of  data  rather  than  single  elements,  which  might  improve 
the  efficiency.  Under  these  conditions,  the  MER  process  is  more  efficient  than 
out-of-core  storage.  Note  that  the  above  procedure  was  done  as  a  test;  in  prac¬ 
tice,  it  is  beneficial  to  use  as  much  of  fast-access  memory  as  is  available,  and 
examples  to  follow  are  based  upon  a  more  efficient  approach. 

In  order  to  judge  the  accuracy  of  the  particular  discretization  employed  for 
this  problem,  data  from  different  levels  of  discretization  are  compared  to  the 
exact  solution  for  a  circular  cylinder  illuminated  by  a  plane  wave.  Table  4.1 
shows  the  convergence  of  the  monostatic  radar  cross  section  as  the  order  of 
discretization  is  refined.  (A  definition  of  radar  cross  section  for  this  problem 
is  available  in  [57].)  Tables  4.2  and  4.3  show  the  current  density  in  the  center 
of  the  shadow  region  and  at  the  specular  point  on  the  cylinder.  The  numerical 
solutions  appear  to  be  converging  to  the  exact  values,  (in  this  and  other 
examples  to  follow,  in  order  to  show  numerical  convergence,  it  is  necessary  to 
replace  the  look-up  tables  used  for  the  Bessel  functions  for  MER  implementation 
with  more  accurate  values.  The  look-up  tables  for  the  Hankel  function  of  order 
zero  exhibit  a  maximum  error  of  1  percent.)  The  numerical  convergence  observed 
in  Tables  4,1  -  4.3  suggests  that  the  various  approximations  used  within  the 
moment  method  formulation  are  acceptable. 

4.4.  TM-wave  Scattering  by  Dielectric  Cylinders 

The  problem  of  TM-wave  scattering  by  lossy,  inhomogeneous  dielectric  cylin¬ 
ders  was  formulated  by  Richmond  [79]  using  pulse  basis  functions  and  point¬ 
matching.  The  MER  implementation  of  this  problem  is  discussed  by  Sultan  and 
Mittra  [13],  [58],  although  no  execution  time  comparisons  are  provided.  Since  the 
off-diagonal  matrix  elements  are  identical  in  form  to  Equation  (4.1),  an  MER 
implementation  yields  relative  efficiencies  identical  to  those  of  the  perfectly 


TABLE  4.2 


CURRENT  DENSITY  AT  CENTER  OF  SHADOW  REGION  INDUCED  BY  PLANE  TM 
WAVE  ON  CIRCULAR  PEC  CYLINDER  WITH  ONE  WAVELENGTH  CIRCUMFERENCE 


N 

\3Z\'K\ 

4 

0.001109 

179.22 

8 

0.000845 

154.16 

16 

0.000784 

152.36 

32 

0.000773 

152.68 

64 

0.000766 

153.00 

128 

0.000763 

153.17 

Exact 

0.000760 

153.35 

TABLE  4.3 


CURRENT  DENSITY  AT  THE  SPECULAR  POINT  INDUCED  BY  PLANE  TM 
WAVE  ON  CIRCULAR  PEC  CYLINDER  WITH  ONE  WAVELENGTH  CIRCUMFERENCE 


N 

Msl/M 

Zli 
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0.006273 

34.557 

8 

0.006391 

41.261 

16 

0.006302 

41.134 

32 

0.006271 

40.792 

64 

0.006254 

40.567 

128 

0.006245 

40.451 

Exact 

0.006237 
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conducting  cylinder  case  discussed  in  Section  4.3.  Thus,  Figure  4.1  suffices  to 
describe  the  dielectric  cylinder  problem  also. 

To  evaluate  the  accuracy  of  the  moment  method  formulation,  numerical  results 
based  on  several  different  models  are  compared  to  the  exact  solution  for  a  homo¬ 
geneous,  circular  cylinder.  The  models  of  the  cylinder  approximate  the  circular 
cross-section  by  a  superposition  of  square  cells,  with  the  areas  of  the  model 
cross-section  normalized  to  the  area  of  the  desired  cylinder.  Table  4.4  shows 
the  monostatic  radar  cross  section  obtained  from  this  process,  and  Table  4.5 
shows  values  of  the  electric  field  at  the  center  of  the  cylinder.  These  results 
suggest  the  validity  of  the  pulse  basis  and  point  matching  formulation,  although 
they  do  not  necessarily  indicate  the  level  of  accuracy  that  might  be  obtained 
for  a  more  complicated  geometry,  such  a.  one  comprised  of  inhomogeneous 
mater ial . 

4.5.  TM-wave  Scattering  by  Cylinders  Modeled  by  a  Combination  of 
Dielectric  and  Conducting  Materials 

The  problem  of  TM-wave  scattering  from  an  infinite  cylinder  containing  both 
dielectric  and  perfectly  conducting  (PEC)  materials  can  be  attacked  by  combining 
the  techniques  discussed  in  the  previous  two  sections.  The  electric  field 
integral  equation  may  be  written  as 
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where  the  first  integral  is  to  be  taken  around  the  outer  surface  of  the  perfect 
conductor  and  the  second  integral  throughout  the  cross-sectional  volume  of  the 
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TABLE  4.4. 

MONOSTATIC  RCS  OBTAINED  FOR  A  HOMOGENEOUS,  CIRCULAR  DIELECTRIC  CYLINDER 
WITH  er  =  10,  AND  CIRCUMFERENCE  OF  0.5137  XQ 


N 

RCS  (dB  Xq) 

21 

-1 .8484 

61 

-1 .8469 

101 

-1 .8442 

exact 

-1.8426 
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dielectric  material.  Note  that  E^  vanishes  for  points  (x,y)  on  the  surface  of 
the  perfect  conductor.  Figure  4.2  illustrates  the  type  of  geometry  under  con¬ 
sideration.  In  Equation  (4.2)  and  throughout  this  report,  k  is  the  free-space 
wavenumber,  n  is  the  intrinsic  impedance  of  free  space,  and  z^  is  the  complex 
relative  permittivity  of  a  dielectric  material. 

A  given  cylinder  can  be  modeled  by  a  superposition  of  N^  perfectly  conducting 
strips  and  homogeneous  dielectric  cylinder  cells.  If  the  current  density  on 
each  strip  and  fields  throughout  each  cell  are  assumed  constant,  and  Equation 
(4.2)  is  enforced  by  point  matching  in  the  center  of  each  strip  and  each  cell, 
the  resulting  system  can  be  written 
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In  Equation  (4.3),  Gdd  denotes  an  N^xN^  matrix,  QPP  denotes  an  N^xNp  matrix,  gdp 

d  "o  es  an  N.xN  matrix,  and  Gpd  denotes  an  N  xN.  matrix.  Ud  and  Up  denote 
a  p  -  p  a  -  - 

arrays  containing  the  coefficients  of  the  pulse  basis  functions  used  to  repre¬ 
sent  the  field  and  current  densities,  which  are  the  unknowns  to  be  determined. 

El  denotes  the  sampled  values  of  the  incident  electric  field. 

If  we  number  the  dielectric  cells  from  1  to  N^  and  the  conducting  strips  from 
Nj+l  to  N^Npi  the  elements  of  the  NxN  matrix,  where  N=Nj+Np,  are  given  by 
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In  Equation  (4.4),  indices  m  and  n  range  from  1  to  .  In  Equation  (4.5),  m  and 
n  range  from  N^+l  to  N^+N^.  G^P  may  be  found  from  Equation  (4.5)  by  treating  the 
index  n  as  ranging  from  N^+l  to  N^+Np,  and  index  m  as  ranging  from  1  to  N^. 
Similarly,  Gpd  may  be  found  from  Equation  (4.4)  with  n  between  1  and  and  m 
between  N^+l  and  N^+Np.  Using  the  approximations  of  Richmond  [79],  the  integral 
in  Equation  (4.4)  can  be  reduced  so  that 
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where  a^  is  the  radius  of  a  circular  cylinder  with  cross-sectional  area  equal  to 

that  of  cell  n,  and  p  is  the  distance  from  the  center  of  cell  n  to  the  center 

mn 

of  cell  m.  Using  a  simple  approximation  suggested  by  Harrington  [74],  the 
integral  of  Equation  (4.5)  can  be  simplified  to  yield 
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where  w  is  the  width  of  the  n-th  strip, 
n 

It  is  not  surprising  that  the  off-diagonal  matrix  elements  are  of  the  form  of 
Equation  (4.1),  since  the  method  is  a  combination  of  two  approaches  which  indi¬ 
vidually  satisfy  this  condition.  The  accuracy  of  these  two  procedures  has  been 
discussed  in  Sections  4.3  and  4.4,  and  additional  convergence  tests  are  not  pro¬ 
vided.  Figures  4.3  and  4.4  illustrate  the  accuracy  of  the  combined  methods  for 
the  surface  current  density  and  internal  fields  found  for  a  circular  PEC 
'•ylinder  coated  with  a  uniform  dielectric  layer. 
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Figure  4.3.  Comparison  of  exact  and  numerical  results  for  the  current  density 
induced  on  a  PEC  cylinder  with  a  dielectric  coating. 
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The  relative  efficiency  of  the  MER  implementation  of  the  above  technique  is 
illustrated  by  Figure  4.5,  which  shows  execution  times  for  both  the  MER  imple¬ 
mentation  and  one  where  matrix  elements  are  transferred  from  a  disk  storage 
device.  In  both  cases  a  125x125  portion  of  the  matrix  was  stored  in  fast-access 
memory  on  the  CDC  CYBER  175,  and  the  remaining  matrix  elements  either  generated 
or  read  from  disk  as  needed.  The  system  was  solved  by  the  iterative  CGM 
algorithm  of  Chapter  2.  Note  that  the  data  shown  represent  the  total  execution 
time,  which  is  proportional  to  the  number  of  iterations  required  for  the  itera¬ 
tive  algorithm  to  converge.  Identical  examples  were  used  for  both  the  MER  and 
the  off-core  storage  data,  and  should  be  a  valid  indication  of  the  relative 
efficiency.  The  total  execution  time  for  other  examples  may  differ  appreciably 
from  these  data  if  a  different  number  of  iterations  are  necessary  to  produce  a 
solution.  From  the  data,  it  is  clear  that  the  MER  approach  is  the  more  efficient 
procedure . 

4.6.  TE-wave  Scattering  by  Cylinders  Modeled  by  a  Combination  of 
Dielectric  and  Conducting  Materials 

A  numerical  solution  to  the  TE-polarization  counterpart  of  the  composite 
dielectric  -  conducting  cylinder  problem  of  Section  4.5  can  be  based  on  a 
variety  of  mathematical  formulations.  Anticipating  an  MER  implementation,  the 
matrix  elements  are  to  be  kept  as  simple  as  possible.  One  of  the  approaches 
involves  the  magnetic  field  integral  equation  (MFIE)  for  the  perfect  conducting 
material  combined  with  the  electric  field  integral  equation  (EFIE)  for  the 
dielectric  material.  An  EFIE  could  be  used  to  treat  the  conducting  material,  but 
would  require  more  complicated  matrix  elements.  Similarly,  the  MFIE  has  been 
applied  to  the  dielectric  material  for  this  polarization,  but  simple  discretiza¬ 
tions  appear  to  yield  more  accurate  results  when  used  with  the  EFIE  under  simi¬ 
lar  cond  i  t ions  [82 ) . 
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The  coupled  system  can  be  written 
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where  A  is  the  magnetic  vector  potential  defined  for  the  conducting  material  as 


A  ( x,y)  =  /  J(l')  )-r  h!:2)  (kR)  dS,’ 
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and  for  the  dielectric  material  as 


Ad(x,y)  =|^  //  (er  -  1)  Et0t  Hg2)(kR)  dx’dy’  (4.10) 

where 

R  =  /(x  -  x’)2  ~  (y  -  y’)2  (4.11) 

Note  that  the  EFIE  of  Equation  (4.8)  is  to  be  enforced  throughout  the  dielectric 
material,  while  the  MFIE  is  to  be  enforced  over  the  surface  of  the  conducting 
material . 

If  the  geometry  to  be  analyzed  is  modeled  exactly  as  in  Section  4.5,  and 
expansion  functions  identical  to  those  of  Section  4.5  are  used  for  the  unknown 
quantities,  Equation  (4.8)  can  be  enforced  at  the  center  of  each  strip  and  cell 
to  yield  a  system  of  the  form 
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If  the  model  contains  dielectric  cells  and  conducting  strips,  there  will 

be  a  total  of  (211^+11^)  unknowns.  We  number  the  x  component  of  the  fields  in  the 

dielectric  cells  from  1  to  ,  the  y  component  of  the  fields  in  the  dielectric 

cells  from  N^  +  '.  to  2N^ ,  and  the  current  density  on  the  conducting  strips  from 

2N,+i  to  2N.+N  .  The  matrix  elements  G,,  have  been  evaluated  approximately  by 
ddp  dd 

Richmond  [80] ,  and  are  given  as 
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for  the  off-diagonal  elements,  where 


Y  =  (er  -  1) 
n  r  n 


j  itka  J .  ( ka  ) 
J  n  1  n 


(4.16) 


X  =  k(x  -  x  ) 
m  n 


(4.17) 


106 


V  =  k(ym  -  yn) 
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Hq(R)  =  Hq2)(R) 
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In  each  of  the  above  expressions  the  indices  refer  to  the  proper  range  of  para¬ 
meters  for  that  sub-matrix.  The  diagonal  elements  are  given  by 


(2) 

,  .  .  .  jitka  H.  (ka  ) 

dd  _  dd  r J  ml  m  ,  . 

G  _  —  G  _  -  1  +  ( £r  —  1 )  -  ■  -  ■]■  -  ♦  1 

xx  mm  yy  mm  rm  L  4 


(4.22) 


d  d  dd 

and  the  diagonals  of  the  sub-matrices  G  and  G  vanish.  The  elements  in  the 

=xy  -yx 

and  g  ^  sub-matrices  can  be  found  by  a  simple  extension  of  Richmond's 
formulas  [80]  to  give 


A  wka  Y 

GxPmn  =  -(£rn-  ^  IT  Jl(kan)  VR)  R 


(4.23) 


,pd 


nka 


C*"*  =  (er  -  D  TT"  J,(kan)  H.(R)  i 

y  mn  cn  2  n  1  n  1  R 


(4.24) 


Note  that  the  m-index  in  Equation  (4.23)  is  to  range  from  (2N.+1)  to  (2N.+N  ) 

d  dp 


and 


the  n-index  from  1  to  .  The  n-index  in  Equation  (4.24)  ranges  from  N^  +  l  to 


2Nd. 

In  order  to  evaluate  the  matrix  elements  for  £PP 
an  approach  used  by  Harrington  [83].  A  single  strip 


G  dp.  and  G  dp,  we  consider 
=x  =y  ’ 

is  depicted  in  Figure  4.6. 
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Note  that  we  need  to  introduce  an  orientation  parameter  71,  which  was  not  needed 
for  the  treatment  of  the  TM  polarization  in  Section  4.5. 

If  the  strips  in  the  model  are  fiat,  the  diagonal  elements  of  Gpp  are 

Gpp  =  -  (4.25) 

mm  2 


The  off-diagonal  elements  are  given  by  Harrington  [83]  as 


Gpp  = 
mn 


kw 


H.UOfi  cos  71 


X 

-  sin 


n 

n 


1 

J 


(4.26) 


Equation  (4.26)  is  actually  an  approximation  of  the  expression 


kw 


GPP  =  i  J  .  H,  (R)  {—  cos  a  -  -  sin  7!  |  dL 
mn  4  J  kw  1  n  e-  n' 


(4.27) 


where 


Y  *  ky  -  ky  -  L  sin  71 
a  7  n  n 


(4.28) 


X  -  kx  -  kx 
m  n 


L  cos  Q 

n 


(4.29) 


and 


(4.30) 


Similarly,  we  find 
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Gdp 
x  mn 


Gdp 
y  mn 


Gdp 
x  mn 


Gdp 
y  mn 


kw 


nr  r  ~  ~  Ho(R)  R  -  2H. (R) 

4  J  ^  tsin  X  Y 


4  J  kw  1  n 
n 


cos  anlH0(R)  - 
R 


?\2 


~  X2  -  Y2 

+  VR)  iL~qL-]}  dL 

n J  * 


7 


(4.31) 


kw 


-  J2  r 
4  J  kw 


r  XY  —  ^ 

{cos  an  -3J  [H0(R)  R  -  2H  (R) ] 

R 


2  2 

■sin  *\[H0(R)[4  +  H  (R)  7  ~  X  jl  dL 

\R/  r3 


(4.32) 


nkwn  .  H„(R)  R  -  2H,(R) 

-r —  sin  t!  XY  - : _ 

4  n 


R" 


cos  an[H0(R)  ~ 


x2  -  Y2 
*  H  (R)  - - JL 

RJ 


(4.33) 


nkwn  ,  H  (R)  R  -  2H,(R) 

-  cos  SI  XY  — - -  1 

4  n 


R 


3 


Sin  WR)(f|  *  Hi(R)  —  ~/2. 

R 


(4.34) 


no 


The  right-hand  side  of  Equation  (4.12)  represents  the  sampled  incident  field 
at  the  centers  of  the  cells  and  strips.  The  unknown  column  array  contains  coef¬ 
ficients  of  the  basis  functions  (pulses  in  this  case)  used  to  represent  the 

fields  in  the  dielectric  and  the  current  density  on  the  conductor. 

Due  to  the  relative  simplicity  of  the  formulas  given  in  Equations  (4.26), 
(4.33)  and  (4.34),  as  compared  to  the  exact  expressions  which  would  normally 
require  some  form  of  numerical  integration  for  evaluation,  they  are  the  initial 
choices  for  MER  implementation.  If  based  upon  these  approximate  expressions,  the 
MER  approach  is  more  efficient  than  transferring  needed  data  from  disk.  As  an 
example,  a  circular  perfectly  conducting  cylinder  coated  by  a  homogeneous 
dielectric  layer  was  described  by  a  150x150  matrix  equation.  If  a  125x125  por¬ 
tion  of  this  equation  is  stored  in  fast-access  memory  on  the  CDC  CYBER  175  com¬ 
puter,  and  the  CGM  algorithm  of  Chapter  2  is  used  to  solve  the  system,  the  MER 

approach  is  approximately  three  times  faster  per  iteration  than  reading  the 
necessary  matrix  elements  from  disk. 

The  time  necessary  to  transfer  needed  data  from  disk  appears  to  be  relatively 
repeatable,  and  for  the  preceding  example  about  4.3  seconds  per  iteration  were 
necessary  to  perform  the  transfers  and  other  computations.  The  MER  implemen¬ 
tation  required  approximately  1.3  seconds  per  iteration.  For  the  TM  polarization 
as  discussed  in  Section  4,5,  an  example  involving  a  150x150  matrix  equation, 
with  a  125x125  part  stored  in  fast  access  memory,  requires  about  0.75  second 
per  iteration  for  MER  implementation.  Clearly,  the  efficiency  of  the  MER 
approach  is  affected  by  the  additional  complexity  of  the  matrix  elements  for  the 
TE  polarization. 

Furthermore,  experimentation  indicates  that  the  approximate  formulas  are  not 
very  accurate  when  the  model  under  consideration  involves  close  spacings  between 
conducting  strips  and  dielectric  cells.  Table  4.6  illustrates  the  accuracy  in 
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TABLE  4.6. 

MAGNITUDE  OF  SURFACE  CURRENT  DENSITY  INDUCED  ON  PEC  CIRCULAR  CYLINDER  BY  PLANE 
TE-POL  WAVE,  WHEN  THE  CYLINDER  IS  COATED  WITH  A  HOMOGENEOUS  LAYER  OF  DIELECTRIC. 

r in  =  °-05  V  rout  =  °-0875  V  £r  =  3'j°’  Np  =  12>  Nd  =  12 

PEC  CELL  DENSITY:  38  strips/Xn 

DIELECTRIC  CELL  DENSITY:  247  cells/X^ 

Numerical  values  produced  using  the  approximate  Equations  (4.26),  (4.33),  and 
(4.34)  compared  to  the  values  obtained  by  a  numerical  integration  of  Equations 
(4.27),  (4.31),  and  (4.32).  The  exact  eigenfunction  solution  is  shown  for  com¬ 
parison  . 


* 

| japp  elem  jy jginc  | 

wm 

o3 

1.776 

1.114 

1.106 

30 

1.719 

1  .075 

1.064 

60 

1.629 

1 .018 

1  .000 

90 

1.670 

1  .059 

1.040 

120 

1 .874 

1.209 

1.198 

150 

2.094 

1  .  363 

1.358 

180 

2.184 

1 .425 

1 .422 
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the  surface  current  density  (magnitude  only)  for  a  circular  PEC  cylinder  coated 
with  a  uniform,  homogeneous  layer  of  dielectric.  In  this  case,  large  errors 
occur  in  the  solution  due  to  the  use  of  the  approximate  formulas  of  Equations 
(4.26),  (4.33),  and  (4.34).  The  accuracy  of  the  approximations  can  be  investi¬ 
gated  based  on  direct  numerical  calculations,  and  for  illustration  Table  4.7 
shows  the  values  of  Equations  (4.33)  and  (4.31)  for  a  range  of  observation 
coordinates.  It  appears  that  the  approximate  expression  is  reasonably  accurate 
except  when  the  observation  point  is  located  within  a  circle  of  radius  0.2  X^, 
centered  at  the  source  strip.  For  close  spacings,  other  approximat ions  or  some 
form  of  numerical  integration  should  be  used  to  evaluate  the  expressions  given 
in  Equations  (4.31)  and  (4.32),  to  ensure  accurate  solutions.  If  implemented 
carefully,  the  more  accurate  expressions  should  not  significantly  interfere  with 
the  efficiency  of  the  MER  approach. 

4.7.  Summary 

This  chapter  has  explored  the  iterative  MER  approach  for  the  solution  of 
electromagnetic  scattering  problems.  The  MER  approach  requires  a  portion  of  the 
system  matrix  to  be  recomputed  whenever  needed  by  the  iterative  algorithm  in 
use,  and  has  the  advantage  that  the  matrix  does  not  need  specific  symmetries  for 
effective  iterative  implementation.  Because  the  procedure  is  not  based  upon  sym¬ 
metries  in  the  matrix  equation,  it  is  not  as  efficient  as  the  other  iterative 
approaches  in  use.  Experimentation  indicates  that  the  MER  approach  can  be  more 
efficient  than  transferring  needed  matrix  elements  from  disk  whenever  needed  in 
the  solution  process,  and  thus  may  be  a  favorable  alternative  for  the  solution 
of  large  matrix  equations.  However,  the  efficiency  of  the  MER  process  depends 
on  the  simplicity  of  the  matrix  elements,  and  will  vary  with  the  problem  under 
cons  iderat ion . 
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TABLE  4.7. 

COMPARISON  OF  NUMERICAL  VALUES  OF  EQUATION  (4.31)  AND  EQUATION  (4.33) 
FOR  A  STRIP  WITH  ORIENTATION  0=0,  wn  =  0 . 1  XQ,xn  =  0,  yn  =  0. 


Eq. 


(4.31) 


Eq . 


(4.33) 


%  error 
in  Eq.  (4.33) 


0.1 

0.07 

0 


0 

0.07 

0.1 


28.03  +  j 147 -8 
26.67  -  j 1 5 .92 
25.21  -  j78 .46 


28.13  +  j 1 14 . 5 
26.78  +  j8 . 36 
25.32  -  j98 . 23 


22% 

78% 

24% 


0.2 

0.14 

0 


0 

0.14 

0.2 


24.03  +  j  29 . 16 
19.12  -  j 9 . 2 52 
13.83  -  j41 . 13 


24.10  +  J27.31 
19.19  -  j 7 . 5 1 0 
13.91  -  j42 .77 


5% 

8% 

4% 


0.3 

0.21 

0 


0 

0.21 

0.3 


18.21  ♦  j5.920 
8.884  -  j 14 . 87 
-1.087  -  j 34 .22 


18.25  +  j  5 . 44 1 
8.923  -  j 14 . 51 
-1.052  -  j34 . 67 


2.5% 

2.1% 

1.3% 


0.5 

0.35 

0 


0 

0.35 

0.5 


5.415  -  j 6 .614 
-8.707  -  jlO.10 
-23.28  -  j 1 2 . 53 


5.358  -  j 6 . 7 5 1 
-8.725  -  j 10 .05 
-23.33  -  j 1 2 . 67 


2% 

0.5% 

0.6% 


1 

0.7 


0 

0.7 


-1.987  +  j  2 . 20 1 


1.7% 
0.7% 
0 


0 


1 


6.035  ♦  j 7 . 199 
15.04  +  jll.25 


-1.999  «•  J2.250 
6.100  +  j 7 . 210 
15.02  +  jll  .30 


.3% 
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Specific  examples  of  cylinders  containing  both  perfectly  conducting  and 
dielectric  materials  were  studied  for  the  MER  approach.  The  TM  polarization 
involves  relatively  simple  off-diagonal  matrix  elements,  and  appears  to  be  a 
good  candidate  for  the  MER  approach.  The  TE  polarization  requires  more  compli¬ 
cated  matrix  elements,  and  the  numerical  efficiency  is  not  as  good  as  that  of 
the  TM  polarization.  Furthermore,  for  the  specific  formulation  used  with  the  TE 
polarization,  care  must  be  taken  when  evaluating  the  matrix  elements  since 
simple  approximations  are  poor  if  the  model  requires  close  spacings  between 
dielectric  and  conducting  cells.  Although  more  accurate  evaluation  of  the  matrix 
elements  (when  required  by  close  spacings)  reduces  the  efficiency  of  the  MER 
procedure,  the  approach  is  still  a  viable  alternative  to  storing  the  needed 
matrix  elements  on  disk. 

Since  the  efficiency  of  the  MER  approach  depends  on  the  complexity  of  the 
matrix  elements,  candidates  for  MER  implementation  must  be  evaluated  on  an 
individual  basis.  Clearly,  the  complexity  of  some  of  the  state-of-the-art 
techniques  for  the  numerical  solution  of  scattering  problems  might  preclude 
their  use  in  the  MER  procedure.  However,  simple  approaches  which  offer 
reasonable  accuracy  should  be  good  candidates  for  MER  solution.  The  MER 
approach  should  be  considered  an  alternative  for  the  treatment  of  large  scat- 
terers  that  do  not  possess  the  necessary  symmetries  to  be  treated  efficiently  by 
other  techniques. 
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5.  A  COMPARISON  OF  TWO  PROCEDURES  FOR  THE  DISCRETIZATION 
OF  CONVOLUTIONAL  INTEGRAL  EQUATIONS 

5.1.  Introduction 

Electromagnetic  scattering  problems  can  often  be  described  by  convolutional 
integral  equations.  These  are  equations  having  the  form 
a 

E(x)  =  /  J(x')  K(x  -  x')  dx'  a  <  x  <  b  (5.1) 

b 

where  J(x)  is  the  unknown  function  to  be  determined.  An  approximate  solution  for 
J(x)  can  be  obtained  by  replacing  Equation  (5.1)  by  a  finite-dimensional 
discrete  system.  If  the  convolutional  form  is  preserved  by  the  discretization 
process,  the  discrete  system  can  be  expressed 

em  =  in  gm-n  m  =  1  >2 . N  (5‘2) 

n=l 


We  adopt  the  notation  of  using  lower-case  letters  to  denote  sequences  and  upper¬ 
case  to  denote  functions.  The  discrete  system  of  Equation  (5.2)  may  be  written 
as  3n  N-th  order  matrix  equation 
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Note  that  the  NxN  matrix  is  Toeplitz,  so  that  each  element  is  a  repetition  of 
the  first  row  or  column  according  to  Equation  (5.4).  Thus,  there  is  a  con¬ 
siderable  amount  of  redundancy  present  in  the  system.  Because  of  this  redun¬ 
dancy,  a  Toeplitz  equation  is  one  example  of  the  type  of  system  for  which 
iterative  solutions  are  well  suited.  It  is  not  a  particularly  good  example, 
however,  because  direct  methods  for  solving  Toeplitz  systems  are  also  available 
and  may  be  more  efficient  than  iterative  algorithms  [54],  [55],  [84]  -  [86].  In 
practice,  convolutional  integral  equations  representing  electromagnetic  scat¬ 
tering  problems  can  often  be  converted  to  discrete  systems  containing  Toeplitz 
or  almost-Toepl itz  symmetries.  An  example  of  an  almost-Toepl itz  system  is  a 
matrix  with  Toeplitz  symmetries  everywhere  except  along  the  main  diagonal,  as 
might  arise  in  connection  with  integral  equations  describing  an  inhomogeneous 
geometry.  These  are  not  usually  well-suited  for  direct  solution,  but  are  easily 
treated  efficiently  using  iterative  algorithms.  Additional  examples  of  the 
iterative  solution  of  almost-Toepl itz  systems  are  presented  in  Chapters  6  and  7. 

Two  types  of  discretization  procedure  have  been  used  in  the  recent  past  to 
convert  convolutional  integral  equations  to  matrix  equations  with  discrete- 
convolutional  symmetries.  Because  of  the  utility  of  these  procedures  in  connec¬ 
tion  with  iterative  solution  algorithms,  a  firm  understanding  of  their 
implementation  is  of  central  importance.  The  first  technique  to  be  considered 
is  the  discrete-convolutional  method  of  moments  (DCMoM),  which  is  a  special  form 
of  the  general  moment  method  procedure  [12],  The  second  is  the  spectral-domain 
fast-Fourier  transform  (SDFFT)  method,  which  is  the  name  we  choose  to  denote  the 
discretization  used  in  connection  with  the  original  "spectral-iterative"  tech¬ 
nique  (SIT)  [111,  [44]  -  [46].  In  contrast  with  the  DCMoM,  which  is  well 
understood  after  years  of  research  on  the  method  of  moments,  the  SDFFT  discreti¬ 
zation  has  not  received  much  exposition  to  date  in  the  literature.  This  chapter 
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presents  the  SDFFT  procedure  in  a  general  form,  so  that  an  equivalence  can  be 
made  between  it  and  the  DCMoM  approach.  It  will  be  shown  that  the  two  procedures 
can  be  identical  in  principle,  although  in  practice  the  SDFFT  can  only  be  con¬ 
sidered  an  approximation  to  the  DCMoM  approach.  Guidelines  for  using  the  SDFFT 
discretization  are  given. 

5.2.  The  Discrete-convolutional  Method  of  Moments  (DCMoM)  Procedure 

Consider  a  convolutional  integral  equation  of  the  form  of  Equation  (5.1).  E 
and  K  are  known  over  the  interval  of  interest  and  J  is  an  unknown  function  to  be 
determined.  Equation  (5.1)  can  be  used  to  describe  scattering  from  a  strip  or 
wire  of  constant  curvature,  and  :  representative  of  a  variety  of  other  electro¬ 
magnetic  scattering  problem*'.  A  discretization  of  Equation  (5.1)  according  to 
the  moment-method  procedure  requires  that  J  be  replaced  by  a  finite  expansion  of 
the  form 


N 

J( x)  -  i  j  B  ( x) 
•  n  n 


(5.5) 


where  the  {B  (x)}  are  known  basis  functions  and  the  j  unknown  coefficients.  If 
n  n 

the  expansion  is  substituted  into  Equation  (5.1)  and  the  resulting  equation  is 
made  orthogonal  to  N  independent  testing  functions  {Tm(x)},  the  result  is  a 
matrix  equation  of  the  form 


e 

m 


N 

L 

n=l 


^  n  ^m , n 


m  =  1 ,2,...  ,N 


(5.6) 


where 


b 

e  =  /  T  (x)  E( x )  dx 
m  m 


(5.7) 
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b  b 

g  =  /  T  (x)  /  B  (x* )  K(x  -  x')  dx'dx  (5.8) 

m ,  n  1  in  J  n 
’  a  a 

In  the  general  case,  gmn  represents  a  fully-populated  matrix  whose  NxN  entries 
satisfy  no  symmetry  or  redundancy  condition. 

If  the  choice  of  basis  and  testing  functions  is  restricted  to  the  form 


Bn(x)  =  B(x  -  xn)  (5.9) 

T  (x)  =  T(x  -  x  )  (5.10) 

m  m 

where 


x  =  xA  +  nAx 
n  0 


n  =  1,2,  ...,N 


(5.11) 


and  if  the  basis  and  testing  functions  do  not  overlap  the  endpoints  of  the 
interval  (a,b),  the  discrete  system  described  in  Equation  (5.6)  can  be  written 
as 


N 

e  =  v-  j  g 
ra  u,  n°m-n 
n=l 


n  »  1,2 . N 


(5.12) 


The  "g"  appearing  in  Equation  (5.12)  is  completely  described  by  only  (2N-1) 
entries  of  the  NxN  system,  and  often  symmetry  considerations  reduce  the  number 
of  independent  entries  to  N.  This  system  is  exactly  the  Toeplitz  form  described 
in  Equations  (5.2)  -  (5.4).  The  moment  method  application  embodied  in  Equations 
(5.9)  -  (5.12)  is  denoted  the  discrete-convolutional  method  of  moments  (DCMoM), 
because  the  summation  appearing  in  Equation  (5.12)  is  a  discrete  convolution. 
Examples  of  the  DCMOM  are  given  in  Chapters  6  and  7,  and  elsewhere  in  the 
literature  [12],  [47]  -  [51]. 
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5.3.  The  Spectral-domain  Fast-Fourier  Transform  (SDFFT)  Method 

The  DCMoM  procedure  is  one  type  of  discretization  which  is  readily  compatible 
with  iterative  solution  algorithms.  A  different  approach,  applicable  to  con¬ 
volutional  integral  equations,  requires  the  Green's  function  K  appearing  in 
Equation  (5.1)  to  be  sampled  in  the  Fourier  transform  or  "spectral"  domain.  In  a 
multidimensional  case,  the  Fourier  transform  is  to  be  taken  with  respect  to  one 
or  more  of  the  spatial  variables.  The  motivation  for  this  alternate  approach 
stems  from  the  fact  that  often  the  Fourier  transform  of  K  is  much  simpler  and 
easier  to  compute  numerically  than  K  itself.  For  instance,  problems  involving 
planar,  stratified  media  give  rise  to  spatial  Green's  functions  in  terms  of 
infinite  integrals  [87]  or  infinite  summations  [88].  If  written  in  the  transform 
domain,  these  Green's  functions  usually  become  algebraic  expressions  [89], 
Unfortunately,  in  practice  many  of  the  integral  equations  of  interest  are  only 
valid  over  a  region  of  finite  support,  and  cannot  be  transformed  entirely  into 
the  spectral  domain.  Instead,  the  fast-Fourier  transform  algorithm  is  used  to 
connect  the  two  domains.  Although  in  the  past  the  procedure  was  developed  from 
a  different  perspective  [44],  it  can  ho  considered  as  nothing  more  than  an 
alternate  way  to  produce  the  matrix  "g"  employed  in  Equation  (5.12),  by  applying 
the  inverse  FFT  to  the  sequence  obtained  by  sampling  the  analytical  Fourier 
transform  of  K(x).  We  denote  this  approach  the  spectral-domain  fast-Fourier 
transform  (SDFFT)  discretization. 

Because  the  SDFFT  approach  requires  a  mixture  of  analytical  Fourier  trans¬ 
form  techniques  and  numerical  applications  of  the  FFT  algorithm,  it  is  important 
to  note  the  differences  beteeen  these  two  tools.  A  detailed  discussion  may  be 
found  in  Brigham  [42].  Specifically,  the  FFT  is  equivalent  to  the  Fourier  trans¬ 
form  only  when  the  latter  is  applied  to  functions  which  are  discrete  and 
periodic  in  both  the  original  and  transform  domains.  As  an  aside,  for  this 
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reason  the  SDFFT  procedure  appears  well-suited  for  the  discretization  of 
equations  which  are  already  periodic  in  one  or  more  spatial  variables.  Some 
examples  of  periodic  equations  in  electromagnetics  arise  in  the  analysis  of 
antenna  arrays  [90]  and  frequency  selective  surfaces  [91].  A  discussion  of  the 
SDFFT  applied  to  periodic  problems  is  reserved  for  a  later  section  of  this 
chapter.  The  topic  of  interest  here  is  the  application  of  the  SDFFT  to  non¬ 
periodic  problems,  for  it  is  in  these  cases  that  additional  care  must  be  taken 
to  ensure  the  accuracy  of  the  process. 

In  order  to  study  the  SDFFT  approach  in  detail,  we  will  construct  the  type 
of  function  for  which  the  FFT  and  Fourier  transforms  are  equivalent,  and 
describe  the  method  using  the  analytical  Fourier  transform.  The  Fourier  trans¬ 
form  is  defined 

oo 

F  { H(  x) }  =  H(f)  =  /  II(x)  e'j2,Tfx  ax  (5.13) 

—  OO 

and  the  inverse  transform  as 

30 

F‘lrH(f)}  =  H(x)  =  /  H(f)  ej2,Tfx  df  (5.14) 

—  <30 

We  adopt  the  conventional  practice  of  denoting  convolution  as 

OO 

A(x)  *  B(x)  =  /  A(x')  B(x  -  x')  dx ’  (5.15) 

—  GO 

Since  the  FFT  is  equivalent  to  the  Fourier  transform  of  a  discrete,  periodic 
function,  it  is  necessary  to  convert  continuous  aperiodic  quantities  to  discrete 
periodic  quantities.  Consider  the  functions 

00 

6(x  -  mAx) 


S(x) 


m=-°o 


(5.16) 
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P(x)  =  V  6( x  -  q  AX)  (5.17) 

q  =  -oo 

Multiplication  with  S(x)  is  equivalent  to  sampling  at  intervals  of  Ax;  con¬ 
volution  with  P(x)  produces  a  periodic  function  with  period  AX.  The  Fourier 
transforms  of  the  above  are 

(5.18) 

(5.19) 

(5.20) 

(5.21) 

In  the  transform  domain,  convolution  with  S( f)  produces  a  periodic  function  with 
period  AF;  mul t ipl icat ion  with  P( f)  produces  a  discrete  function  sampled  at 
intervals  Af.  In  practice,  the  periods  and  sampling  intervals  are  related  by  an 
integer  M,  so  that 

AX  =  M  Ax  (5.22) 

AF  =  M  Af  (5.23) 


S( f)  =  AF  V  6( f  -  mAF) 
m =  00 


P(f)  »  Af  l  S(f  -  q Af ) 

q  =  -oo 


where 


AF  = 


Ax 


Af  =ix 


Equations  (5.15)  -  (5.19)  will  be  used  to  convert  functions  to  discrete, 
periodic  sequences  in  order  to  model  the  FFT  algorithm. 

According  to  the  discussion  of  the  SIT  in  the  literature  [11],  the  SDFFT 
procedure  involves  sampling  the  transform  K( f )  in  such  a  manner  as  to  create  a 
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discrete,  periodic  "spectral  Green's  function"  of  the  form 


Gj ( f )  =  S(f)  *  tP(f)  W(f)  K(f)]  (5.24) 

where  K( f )  represents  the  analytical  Fourier  transform  of  K(x),  S(f)  and  P(f) 
are  defined  above,  and  W( f)  is  a  windowing  function  introduced  to  truncate  the 
support  of  K( f ) ,  typically  to  one  period.  This  construction  allows  us  to  imme¬ 
diately  transform  the  discrete  spectral  Green's  function  to  the  spatial  domain, 
and  make  a  comparison  with  the  analogous  quantity  arising  from  the  DCMoM  proce¬ 
dure.  The  discrete  spatial  Green's  function  obtained  from  Equation  (5.24)  is 
given  by 

Gj(x)  =  P(x)  *  (S(x){W(x)  *  K.  ( x )  }  ]  (5.25) 


An  alternate  way  of  denoting  this  is 


g 


(  1) 

H-m 


'l  W(  x)  *  K(  x) 
q  =  -a> 

X 


( £-m) Ax-q AX 


(5.26) 


If  the  DCMoM  process  is  generalized  to  produce  an  infinite-periodic 
sequence  that  coincides  with  the  previous  DCMoM  sequence  "g"  over  the  interval 
of  interest  in  the  spatial  domain,  a  similar  procedure  can  be  used  to  produce 
the  discrete  Green's  function 


G2(x)  =  P(x)  *  (S(x)  U( x) (T( -x)  *  B( x )  *  K( x)  }  ] 


(5.27) 


or,  equivalently, 
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g 


(2) 

£-m 


T(-x) 


B(x)  *  K( x ) 


(5.28) 


x  =  (jJ-m)Ax 

B(x)  is  the  basis  function  introduced  in  Equation  (5.9),  and  T(-x)  is  a  space- 
reversal  of  the  testing  function  T(x)  appearing  in  Equation  (5.10).  The  notable 
difference  between  the  functional  form  of  Equations  (5.25)  and  (5.27)  is  the 
appearance  of  U(x)  in  the  DCMoM  function.  U(x)  is  necessary  to  truncate  the  spa¬ 
tial  kernel  K(x)  to  the  period  in  order  to  avoid  aliasing  errors.  For  instance, 
U(x)  may  be  of  the  form 


U(x) 


x  £  (a,b) 
otherwise 


(5.29) 


Of  course,  the  period  may  be  larger  than  the  interval  (a,b),  and  U(x)  may  vary 
accordingly.  The  aliasing  errors  due  to  the  absence  of  U(x)  are  clearly 
illustrated  by  the  infinite  summation  in  Equation  (5.26). 

In  view  of  the  above  comparison,  it  appears  that  the  SDFFT  process  should 
be  generalized  to  incorporate  a  function  corresponding  to  the  U(x)  used  with  the 
DCMoM  technique.  However,  U(f)  appears  within  a  convolution  in  the  spectral 
domain,  and  the  U( f )  corresponding  to  Equation  (5.29)  is  a  so-called  sine  func¬ 
tion,  with  support  over  the  entire  x-axis.  Because  of  this,  in  general  it  is 
difficult  to  include  the  convolution  with  U(f)  in  a  numerical  implementation. 

There  appear  to  be  two  ways  in  which  the  effects  of  U(x)  could  be  included 
approximately  in  the  SDFFT  procedure.  The  first  is  simply  to  extend  the  period 
to  some  large  interval,  and  approximate  the  transform  U(f)  by  a  Dirac  delta 
function  (which  it  approaches  as  the  period  becomes  sufficiently  large).  This  is 
the  technique  used  in  the  literature  [11],  [13],  [44]  -  [46].  (Note  that  a 
given  function  is  not  altered  after  convolution  with  a  delta  function,  and  thus 
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Equation  (5.25)  suffices  to  describe  the  process.)  An  alternate  approach  is  to 
replace  the  rectangular  truncation  function  of  Equation  (5.29)  with  a  smoother 
function,  in  order  to  obtain  a  U(f)  with  finite  support  (at  least 
approximately).  The  smoother  U(x)  must  be  sufficiently  flat  over  the  spatial 
interval  of  interest,  in  order  to  avoid  distorting  the  desired  spatial  Green's 
function,  yet  yield  a  transform  which  can  be  conveniently  included  in  the  con¬ 
volution  operation  of  the  discrete  spectral  Green's  function 

Gj ( f )  =  S(f)  *  lU(f)  *  P(f)  W(f)  K(f)J  (5.30) 

The  first  approximate  method  for  including  U(x)  has  the  advantage  that  the 
convolution  with  U( f)  disappears  from  Equation  (5.30),  simplifying  the  calcula¬ 
tion  of  G^(f).  However,  the  period  must  be  significantly  larger  than  would 
otherwise  be  the  case,  in  order  to  ensure  that  the  summation  of  Equation  (5.26) 
is  an  adequate  approximation  to  the  desired  result  (the  q=0  term  alone). 

Examples  to  follow  will  attempt  to  determine  the  size  period  necessary  for 
reasonable  accuracy.  It  is  important  to  note  that  the  large  period  is  only 
necessary  for  the  initial  construction  of  G(f);  once  a  suitable  discrete 
spectral  Green's  function  is  obtained  it  can  be  transformed  (via  the  FFT)  to  the 
spatial  domain  and  truncated  to  the  interval  of  interest.  The  second  approach 
has  the  drawback  that  additional  computation  will  be  necessary  to  include  the 
effects  of  the  approximate  U(x),  assuming  such  a  function  can  be  found.  An 
advantage  of  the  second  approach  is  that  a  singularity  often  present  in  K(f)  is 
explicitly  smoothed  by  convolution  with  U(f).  Ray  and  Mittra  have  attributed 
irregular  numerical  results  to  improper  treatment  of  the  singularity  in  K( f) 
(92]. 
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It  is  clear  from  a  comparison  of  Equations  (5.26)  and  (5.28)  that  the  func¬ 
tion  obtained  by  convolving  the  basis  function  B(x)  and  the  space-reversed 
testing  function  T(-x)  in  the  DCMoM  process  plays  an  analogous  role  to  the 
inverse  transform  of  the  windowing  function  W( f )  used  in  the  SDFFT  procedure. 
This  suggests  that  the  SDFFT  is  equivalent  to  the  DCMoM,  if  the  latter  uses 
basis  functions  W(x)  and  Dirac  delta  testing  functions.  This  fact  was  recently 
noted  by  Nvo  and  Harrington  [47].  Previous  presentations  of  the  SDFFT  approach 
in  the  context  of  SIT  have  omitted  the  point  that  the  spatial-domain  basis  func¬ 
tion  is  implicitly  chosen  as  the  inverse  transform  of  the  windowing  function.  Of 
course,  this  equivalency  only  holds  if  the  aliasing  effects  of  the  summation  in 
Equation  (5.26)  are  negligible. 

Because  the  windowing  function  W( f)  appears  as  a  multiplication  with  K( f) , 
the  SDFFT  approach  can  easily  incorporate  a  variety  of  windowing  functions. 

Based  upon  the  comparison  with  the  DCMoM,  it  appears  that  a  primary  con¬ 
sideration  for  the  choice  of  W(  f)  should  be  the  corresponding  spatial  domain 
basis  function  selected  implicitly  in  the  process.  For  instance,  the  choice  of 
a  rectangular  window  for  W( f)  corresponds  to  the  implicit  choice  of  a  sine  func¬ 
tion  for  the  basis  function.  Since  sine  functions  have  unbounded  support,  they 
do  not  appear  to  be  appropriate  approximations  to  subsectional  basis  functions, 
and  will  apparently  have  considerable  support  outside  the  original  domain  of 
interest  (i.e.,  outside  the  original  scatterer).  Typical  subsectional  basis 
functions  such  as  a  piecewise  constant  or  a  triangle  function  thus  correspond  to 
windowing  functions  W(f)  with  unbounded  support,  which  seems  to  suggest  that  the 
"proper"  windowing  function  to  use  with  the  SDFFT  is  one  which  allows  con¬ 
siderable  aliasing  in  the  spectral  domain.  Thus,  the  incorporation  of  a  win¬ 
dowing  function  that  corresponds  to  a  subsectional  basis  function  may  be 


126 


complicated  by  the  need  to  deliberately  overlap  many  periods  of  the  function 
K(f)  when  constructing  G^(f). 

In  practice,  the  selection  of  W(f)  may  be  thought  of  as  a  scheme  to  produce 
a  reasonable  approximation  to  a  subsectional  basis  function  and  a  reasonable 
approximation  to  a  windowing  function  with  finite  support.  For  instance,  Figures 
5.1  and  5.2  show  rectangular  and  exponential  windowing  functions,  and  their 
inverse  transforms.  The  exponential  window  might  be  an  improvement  over  the  rec¬ 
tangular  window  because  the  basis  function  associated  with  the  exponential  win¬ 
dow  has  its  support  reasonably  confined,  yet  is  similar  to  that  produced  by  the 
rectangular  window  over  the  desired  interval.  Furthermore,  the  exponential  win¬ 
dow  only  requires  approximately  two  periods  to  overlap  in  the  calculation  of 

Gj ( f )  . 

By  analogy  with  the  DCMoM  procedure,  it  is  obvious  that  -a  testing  function 
could  be  incorporated  into  the  SDFFT  process,  as  may  be  necessary  if  the  excita¬ 
tion  in  a  given  problem  is  highly  localized.  In  principle,  the  choice  of  W(  f ) 
can  be  made  to  correspond  to  both  a  basis  and  a  testing  function,  and  the  exci¬ 
tation  sequence  "e"  can  be  computed  according  to  Equation  (5.7). 

5.4.  Iterative  Implementation  of  the  DCMoM  and  SDFFT  Systems 

From  the  above  discussion,  it  is  apparent  that  the  DCMoM  system  and  the  SDFFT 
system  can  both  be  represented  by  Equation  (5.12),  with  the  upper  limit  of  the 
summation  limited  to  N.  In  this  case,  N  represents  the  number  of  points  lying  in 
the  original  domain  (a,b).  Note  that  there  are  2N-1  pertinent  values  of  "g" 
appearing  in  Equation  (5.12).  Thus,  both  the  SDFFT  and  DCMoM  systems  can  be 
implemented  in  exactly  the  same  manner.  While  this  implementation  can  be 
accomplished  in  many  ways,  for  instance  by  solving  an  N-th  order  matrix 
equation,  our  interest  here  is  centered  on  the  iterative  algorithms  of  Chapter  2 
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Figure  5.1.  Comparison  of  the  windowing  functions 

n  ifi  <  0.5 


Wj(f)  =  ' 


otherwise 


and 


W  ( f )  =  exp( -{0.7  7 f } ‘ 


W(f) 


12£ 


X 


Figure  5.2.  Comparison  of  the  implicit  basis  functions 
W  ^  (  x )  *  s in( tx) / ( nx) 


and 


W2(x) 
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These  algorithms  require  successive  applications  of  the  operator  (NxN  matrix) 
described  in  Equation  (5.12). 

The  operator  can  be  applied  in  several  ways.  One  approach  is  to  perform  the 
matrix  multiplication  explicitly  whenever  needed,  which  does  not  imply  that  the 
full  NxN  matrix  must  be  stored  in  computer  memory.  A  popular  alternative  is 
to  make  use  of  the  discrete  convolution  theorem  [42],  and  use  the  FFT  algorithm 
to  perform  the  convolution  of  Equation  (5.12),  according  to 

N  2N-1 

V 

L 

n=  1 

=  FFT'Mj  g}  (5.32) 


■^ngm-n  “ . 

n=l 


Jngm-n 


(5.31) 


where  "j"  must  be  zero-padded  according  to 

j  =  0  n=N+  l,N+2,.. . , 2N- 1  (5.33) 

n 

This  approach  typically  requires  2Mlog(M)  operations,  where  M  =  2N-1  and  the 

2 

logarithm  is  with  respect  to  base  2,  as  opposed  to  N  if  the  operator  was 
implemented  as  a  matrix  multiplication.  Thus,  there  are  fewer  computations 
required  using  the  FFT  assisted  approach  for  a  one-dimensional  system  as  long  as 
N  is  greater  than  20.  Two-dimensional  systems  of  the  form 

N  M 

)'  }  i  g  (5.34) 

,  ,  Jnm  6k-n, fc-m 

n=l  m=l 

can  be  treated  in  a  similar  manner,  as  can  three-dimensional  equations.  Note 
that  for  two-dimensional  systems,  the  total  number  of  unknowns  must  be  approxi¬ 
mately  125  before  the  FFT-assisted  approach  exceeds  the  efficiency  of  direct 
matrix  multiplication.  Because  of  the  zero-padding  requirements,  the  FFT 
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approach  requires  larger  array  sizes,  which  also  affects  the  trade-off  between 
the  two  approaches . 

5.5.  Example:  TM-wave  Scattering  by  Strips 

As  an  example  of  the  implementation  of  the  DCMoM  and  SDFFT  procedures,  and 
a  means  to  compare  the  two  to  develop  guidelines  for  the  use  of  the  latter 
approach,  consider  the  problem  of  TM-wave  scattering  by  a  perfectly  conducting 
flat  strip.  The  integral  equation  for  a  one  wavelength  strip  has  the  form 

0.95 

E(x)  =  /  J(x')  K(x  -  x')  dx'  -  0.05  <  x  <  0.95  (5.35) 

-0.05 


The  kernel  in  this  case  is  given  by 


K(x)  =  -jy  Hq  ( 2 tt  ] x  | ) 


and  its  Fourier  transform  is 

r 


1 


K(f)  =/ 


j4rr/l  -  f2 

1 


V 


4^/f2  -  1 


f  <  1 


If  >  1 


(5.36) 


(5.37) 


If  ten  basis  functions  are  used  with  the  DCMoM  procedure,  specifically  piecewise 
constant  or  "pulse"  functions  defined  by 


B(x)  =  \ 


V 


x  C  (-0.05,  0.05) 

otherwise 


(5.38) 


and  if  the  testing  functions  are  Dirac  delta  functions 
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T(x)  =  6 ( x )  (5.39) 

where 

xn  =  -0.1  +  n(0.1)  n  =  1,2,.  ...10  (5.40) 

then  the  numerical  values  of  the  spatial  domain  sequence  representing  "g"  are 
given  in  Table  5.1.  This  sequence  was  obtained  from  Equation  (5.8),  and  is 
equivalent  to  the  values  of  the  discrete  spatial  Green's  function  G^fx) 
described  in  Equations  (5.27)  or  (5.28). 

Consider  an  SDFFT  discretization  of  Equation  (5.35),  assuming  that  the 
function  U( f )  is  approximated  by  a  Dirac  delta  function  as  discussed  in  Section 
5.3.  Initially,  a  rectangular  window  is  used  for  W(f).  The  objective  of  this 
example  is  to  determine  the  value  of  M  from  Equation  (5.22)  that  produces  a 
reasonably  accurate  discrete  spatial  Green's  function  over  the  interval 
(-0.05,0.95),  if  it  is  possible  to  do  so  and  assuming  that  the  values  produced 
by  the  DCMoM  are  interpreted  as  the  "correct"  result.  Note  that  the  data  from 
Table  5.1  are  not  the  numerical  values  which  would  be  produced  by  the  SDFFT  pro¬ 
cedure  even  if  an  infinite  amount  of  zero-padding  was  incorporated  into  the  pro¬ 
cess,  because  B(x)  and  W( f)  are  not  a  transform  pair. 

Tables  5.2,  5.3,  and  5.4  show  values  of  the  sequence  "g"  produced  by  the 
SDFFT  procedure  for  values  of  M  of  99,  255,  and  1023.  These  tables  also  show  the 
relative  difference  between  the  data  of  Table  5.1  and  the  SDFFT  data.  In  spite 
of  the  fact  that  we  do  not  expect  perfect  agreement,  it  appears  that  the  numeri¬ 
cal  values  are  similar  provided  that  M  is  large.  From  a  study  using  a  variety 
of  strip  sizes,  it  appears  that  the  equivalent  spatial  period  necessary  for 
agreement  between  the  DCMoM  sequence  and  the  SDFFT  sequence  must  exceed  25 
wavelengths  in  order  to  obtain  agreement  within  five  percent  in  the  first  few 


TABLE  5.1 


DISCRETE  SPATIAL  DOMAIN  SEQUENCE  PRODUCED  BY  DCMOM  WITH 
PULSE  BASIS  FUNCTIONS  AND  DIRAC  DELTA  TESTING  FUNCTIONS 
FOR  A  STRIP  OF  LENGTH  1  A  WITH  10  CELLS. 


n 

igni 

J gn  (degrees) 

0 

0.0440 

-34.64 

1 

0.0237 

-71.33 

2 

0.0171 

-111.38 

3 

0.0141 

-149.03 

4 

0.0123 

174.06 

5 

0.0110 

137.49 

6 

0.0101 

101.10 

7 

0.0093 

64.81 

8 

0.0087 

28.59 

9 

0.0082 

-7.58 

TABLE  5.2 

DISCRETE  SPATIAL  DOMAIN  SEQUENCE  PRODUCED  BY  SDFFT  USING 
A  RECTANGULAR  WINDOW  W( f)  FOR  M=99 . 


n 

l*J 

J gn  (degrees) 

X  d i f  f  .  as 
compared  to 
Tabic  5.1 

0 

0.0469 

-26.65 

16  X 

1 

0.0226 

-56.75 

25  X 

0.0159 

-110.65 

7  % 

0.0170 

-146.70 

21  % 

0.0172 

-168.59 

54  X 

5 

0.0130 

171.04 

65  % 

6 

0.0078 

130.44 

50  % 

7 

0.0087 

64.03 

7  X 

8 

0.0125 

35.37 

45  X 

9 

0.0130 

20.48 

84  % 

TABLE  5.3 

DISCRETE  SPATIAL  DOMAIN  SEQUENCE  PRODUCED  BY  THE  SDFFT  USING 
A  RECTANGULAR  WINDOW  W(f)  FOR  M=255. 


gn  (degrees) 


%  diff.  as 
compared  to 
Table  5.1 


0.0425 

0.0227 

0.0174 

0.0140 

0.0113 

0.0092 

0.0091 

0.0094 

0.0086 

0.0070 


34 
-72.14 
-114.94 
-146.48 

178.67 

136.67 
95.37 
61.78 
32.21 
-0.69 


3.5  % 
4.2  % 

6.4  % 

4.5  % 
11  % 
16  % 
13  Z 

5.4  % 
8.1  % 
18  % 


TABLE  5.4 

DISCRETE  SPATIAL  SEQUENCE  PRODUCED  BY  THE  SDFFT  USING 
A  RECTANGULAR  WINDOW  W(f)  WITH  M=1023. 


n 

18  1 
n  1 

y/gn  (degrees) 

%  diff.  as 
compared  to 
Table  5.1 

0 

0.0433 

-35.37 

1.9  % 

1 

0.0237 

-72.47 

2.0  % 

2 

0.0177 

-114.32 

6.0  Z 

3 

0.0139 

-148.42 

2.1% 

4 

0.0117 

173.18 

5.3  % 

5 

0.0104 

132.60 

10  % 

6 

0.0101 

96.09 

8.8  % 

7 

0.0097 

62.70 

5.1  % 

8 

0.0085 

29.61 

3.0  % 

9 

0.0075 

-8.95 

— — - 

9.5  % 
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values  of  the  sequence  "g."  The  period  must  exceed  100  wavelengths  ii  the  first 
few  values  of  the  sequence  are  to  agree  within  one  percent. 

In  an  effort  to  make  a  more  systematic  comparison,  consid«-  tue  DCMoM  and 
SDFFT  procedures  applied  to  Equation  (5.35)  when  both  use  the  same  basis  func¬ 
tions.  The  basis  function 


B(  x) 


exp(-  {■ 


X 

0.7  Ax 


2) 


(5.41) 


is  easily  used  within  the  DCMoM  procedure  since  it  can  be  approximated  by  a 
function  of  finite  support  for  numerical  calculations.  Its  Fourier  transform 


B(f)  =  Ax  exp(-  {— : )  (5.42) 

is  also  a  good  approx imat ion  to  a  function  of  finite  support,  and  can  be  used  as 
a  windowing  function  within  the  SDFFT  process.  (These  functions  are  depicted  in 
Figures  5.1  and  5.2.)  Numerical  values  of  the  sequence  "g"  obtained  by  the  SDFFT 
process  with  M  equal  to  255  and  1023  are  presented  in  Tables  5.5  and  5.6,  and 
are  compared  to  the  sequence  produced  by  the  DCMoM  process  using  the  basis  func¬ 
tions  of  Equation  (5.41).  The  relative  accuracy  obtained  with  a  certain  level  of 
2ero-padding  agrees  well  with  the  previous  results  of  Tables  5.3  and  5.4,  and 
supports  the  above  conclusions. 

The  difference  between  the  SDFFT  sequence  and  the  DCMoM  sequence  as 
illustrated  above  appears  to  be  entirely  due  to  the  fictitious  periodic  nature 
of  the  SDFFT  representation  (i.e.,  the  finite  M)  and  not  to  direct  numerical 
difficulties  associated  with  the  singularity  in  K(f).  Of  course,  if  M  is  taken 
to  be  an  integer  multiple  of  AF,  a  sample  point  will  coincide  with  the  singu¬ 
larity  at  |f|*l  and  the  numerical  procedure  will  fail.  As  described  in  Section 
5.3,  it  should  be  possible  to  incorporate  a  truncation  function  into  the  SDFFT 


TABLE  5.5 


SPATIAL  DOMAIN  SEQUENCE  PRODUCED  BY  THE  SDFFT  USING  AN 
EXPONENTIAL  WINDOWING  FUNCTION  W( f )  WITH  M=255.  THE 
RELATIVE  DIFFERENCE  BETWEEN  THE  SDFFT  AND  DCMOM  SEQUENCE 
FOR  THE  SAME  BASIS  FUNCTIONS  IS  SHOWN  (COMPARISONS  NOT 
AVAILABLE  FOR  n=0  AND  n=l) 


n 

lgJ 

y/gn  (degrees) 

dif  ference 
with  DCMoM 

0 

0.0369 

-39.45 

_ 

1 

0.0228 

-68.00 

- 

2 

0.0166 

-111.26 

2.8  % 

3 

0.0136 

-145.30 

4.8  % 

4 

0.0108 

179.86 

13  % 

5 

0.0089 

138.05 

17  % 

6 

0.0087 

95.93 

14  % 

7 

0.0090 

62.61 

4.6  % 

8 

0.0082 

33.70 

8.5  2 

9 

0.0067 

-0.13 

20  % 

TABLE  5.6 

SPATIAL  DOMAIN  SEQUENCE  PRODUCED  BY  THE  SDFFT  USING  AN 
EXPONENTIAL  WINDOWING  FUNCTION  W(f)  WITH  M=1023.  THE 
RELATIVE  DIFFERENCE  BETWEEN  THE  SDFFT  AND  DCMOM  SEQUENCE 
FOR  THE  SAME  BASIS  FUNCTIONS  IS  SHOWN.  (COMPARISONS  NOT 
AVAILABLE  FOR  n=0  AND  n=l). 


n 

!grJ 

/gn  (degrees) 

dif  ference 
with  DCMoM 

0 

0.0378 

-40 .46 

1 

0.0238 

-68.53 

- 

2 

0.0169 

-110.67 

1.9  % 

3 

0.0135 

-147.21 

2.4  Z 

4 

0.0111 

174.35 

1 .2  X 

5 

0.0099 

133.86 

9.9  % 

6 

0.0097 

96.59 

8.7  X 

7 

0.0092 

63.52 

3.4  X 

8 

0.0081 

30.11 

4.1  % 

9 

0.0071 

-8.42 

11  l 

136 


process,  which  will  reduce  the  necessary  zero-padding  and  explicitly  smooth  the 
singularity  in  K(f).  This  may  be  necessary  in  a  multidimensional  problem  in 
order  to  keep  the  necessary  array  sizes  within  the  storage  constraints  of  the 
computer  in  use. 

It  is  interesting  that  the  SDFFT  sequence  based  upon  the  implicit  basis 
funct ion 


,,,  .  sin(  nx/  Ax) 

w(x)  ~  Tiro t . . 


(5.43) 


agrees  well  with  the  sequence  produced  by  the  DCMoM  with  pulse  basis  functions 
(assuming  the  aliasing  effects  are  suppressed).  Since  it  appears  that  most  of 
the  previous  results  obtained  with  the  SDFFT  used  implicit  basis  functions  of 
the  form  of  Equation  (5.43),  this  may  explain  the  reported  success  of  the  proce¬ 
dure  [11],  [13],  [44]  -  [46],  [50]. 

5.6.  Application  of  the  SDFFT  Procedure  to  Periodic  Equations 

Integral  equations  representing  infinite-periodic  geometries  such  as 
idealized  antenna  arrays  [90]  or  frequency  selective  surfaces  [91]  involve  a 
Green's  function  K(x)  that  is  periodic.  If  Equation  (5.1)  represents  a 
periodic  problem,  it  can  be  discretized  using  the  SDFFT  procedure  without  the 
detrimental  effects  introduced  by  the  periodic  nature  of  the  FFT  algorithm.  For 
instance,  since  the  Fourier  transform  of  a  periodic  function  is  discrete, 
Equation  (5.30)  simplifies  to 

Gj ( f )  =  S(f)  *  [W( f)  K( f ) ]  (5.44) 


and  Equation  (5.26)  is  given  by 

H-l  =  w(x)  *  K(x)l 


x=( i-m) Ax 


(5.45) 
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Thus,  the  SDFFT  process  and  the  DCMoM  process  are  equivalent  for  periodic 
equations  if  W(x)  is  chosen  to  be  the  desired  basis  function  (or  bas is-test ing 
function  pair).  Since  the  spatial  Green's  function  for  a  periodic  problem  is 
usually  an  infinite  summation,  it  may  be  easier  to  work  with  the  SDFFT  than  with 
the  DCMoM.  Examples  of  the  SDFFT  applied  to  periodic  equations  are  available  in 
the  literature  [91],  [93]. 

5.7.  Summary 

This  chapter  presents  a  comparison  of  two  procedures  for  the  discretization 
of  convolutional  integral  equations,  the  discrete-convolutional  method  of 
moments  (DCMoM)  and  the  spectral-domain  fast-Fourier  transform  method  (SDFFT). 
Forms  of  the  DCMoM  have  been  used  by  Bojarski  [10],  Nyo  and  Harrington  [12], 
[47],  Borup  and  Gandhi  [48],  [49],  Ray  and  Mittra  [50]  and  Hurst  and  Mittra 
[51].  The  SDFFT  discret izat ion  has  been  used  in  connection  with  the  spectral- 
iterative  technique  (SIT)  developed  by  Ko  and  Mittra  [44],  Kastner  and  Mittra 
[11],  [45],  [46],  Tsao  and  Mittra  [91],  and  Sultan  and  Mittra  [13].  For  multi¬ 
dimensional  problems,  the  SIT  actually  used  a  DCMoM  discretization  in  one 
variable  and  an  SDFFT  discretization  in  the  others  [11]. 

For  problems  involving  planar  stratified  media  or  problems  that  are 
periodic  in  one  or  more  spatial  variable,  the  SDFFT  process  might  be  easier  to 
implement  and  computationally  more  efficient  than  the  DCMoM  because  of  the  ease 
of  working  directly  with  the  Fourier  transform  of  the  Green's  function.  However, 
for  non-periodic  equations,  a  large  amount  of  zero-padding  may  be  required  to 
initially  construct  the  discrete  spectral  Green's  function  from  the  analytical 
transform.  Because  of  this,  it  may  be  difficult  to  ensure  an  arbitrary  degree  of 
accuracy  in  the  discrete  spectral  Green's  function.  Since  the  DCMoM  procedure 
works  directly  in  the  spatial  domain,  it  appears  to  be  preferable  to  the  SDFFT 
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for  problems  where  the  Green's  function  is  readily  available  in  the  spatial 
domain . 

Based  upon  a  comparison  between  the  DCMoM  and  SDFFT  procedures,  it  is 
apparent  that  both  produce  a  matrix  equation  with  Toeplitz  or  almost-Toepl itz 
symmetries  that  can  be  interpreted  as  an  approximation  to  the  integral  equation 
of  interest.  Iterative  algorithms  for  solving  the  SDFFT  equation  can  be  applied 
in  an  identical  manner  to  solve  the  DCMoM  equation,  and  vice  versa.  Both  methods 
should  be  considered  alternative  possibilities,  with  relative  advantages  and 
disadvantages  dependent  on  the  specific  integral  equation  under  consideration. 
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6.  EFFICIENT  ITERATIVE  IMPLEMENTATION  FOR  TE-SCATTERING  BY  DIELECTRIC  CYLINDERS 

6.1.  Introduction 

Preceding  chapters  have  presented  various  aspects  of  iterative  approaches 
for  the  numerical  solution  of  electromagnetic  scattering  problems.  It  was  men¬ 
tioned  that  the  most  efficient  approaches  all  rely  on  discrete-convolutional 
symmetries  in  the  matrix  equation.  Some  indication  of  how  these  symmetries  can 
arise  was  given  in  Chapter  5,  which  discussed  discretization  procedures  such  as 
the  DCMoM.  This  chapter  presents  specific  details  concerning  the  analysis  of 
TE-wave  scattering  by  lossy,  inhomogeneous  dielectric  cylinders.  It  is  shown 
that  a  DCMoM  discretization  of  the  electric  field  integral  equation  can  yield 
a  2x2  block  system,  where  each  of  the  4  blocks  is  itself  a  block-Toeplitz  or 
almost-block-Toepl itz  matrix.  The  low-storage  iterative  implementation  of 
this  problem  requires  the  organization  of  the  matrix  equation  to  best  expLoit 
the  symmetries.  This  chapter  presents  both  the  implementation  of  the  approach 
and  an  evaluation  of  the  accuracy  of  the  overall  numerical  process. 

A  detailed  description  of  the  problem  of  TE-wave  scattering  by  dielectric 
cylinders  is  provided  by  Richmond  [80],  who  used  pulse  basis  functions  and 
point  matching  to  discretize  an  electric  field  integral  equation.  Below,  we 
briefly  review  the  Richmond  formulation  applied  to  a  special  type  of  model.  In 
particular,  if  the  model  under  consideration  is  restricted  to  that  whose  cross- 
section  is  a  lattice  of  evenly-spaced  square  cells,  such  as  depicted  in  Figure 

6.1,  the  resulting  matrix  equation  contains  the  desired  discrete-convolutional 
symmetries.  In  general,  each  cell  of  such  a  model  may  represent  a  region  of  dif 
ferent  permittivity  without  affecting  the  significant  symmetry  features.  Lossy 
regions  may  be  modeled  by  the  use  of  complex-valued  permittivities. 
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6.2.  Formulation  of  the  Matrix  Equation 

Suppose  the  cylindrical  scatterer  depicted  in  Figure  6.1,  which  is  charac¬ 
terized  by  a  complex-valued  relative  permittivity  £  (x,y),  exists  in  the  pre¬ 
sence  of  a  source  J  .  (Note  that  £  =1  outside  the  cylinder.)  If  the 

source  r  3 

scatterer  were  absent,  the  source  J  would  radiate  a  TE  field  E  .  Because 

source 

of  the  presence  of  the  cylinder,  the  actual  field  differs  from  EinC.  The  task  is 
to  find  this  field,  which  we  denote  Et0t:.  These  fields  are  related  at  all  points 
in  space  by  [ 80 ] 

— inc  —tot  curl  curl  A  -  J  ,, 

E  =  E - : -  (6.1) 

Jaeo 

where 


A(x,y)  =  ~r~-  J  J  J  (kR)  dx'dy' 


(6.2) 


■  A 


x-x' )2  +  ( y—  y '  )2 


(6.3) 


J  =  j'je0(F-r  -  1)  Et0t  (6.4) 

The  integration  in  Equation  (6.2)  is  over  the  support  of  the  polarization 
current  density  J,  which  vanishes  outside  the  scatterer. 

If  restricted  to  the  support  of  the  scatterer,  Equation  (6.1)  is  an  integro- 
differential  equation  for  the  unknown  J(x,y).  Once  J  is  determined,  secondary 
quantities  of  interest  can  be  found  from  the  magnetic  vector  potential  A  g:ven 
in  Equation  (6.2),  using  standard  formulas  such  as 


Trtot  _  -Trine  .  .  -r 

H  =  H  +  curl  A 


(6.5) 
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curl  curl  A  -  J 

Etot  =  Elnc  + - - -  (6.6) 

J“eo 

Since  Equation  (6.1)  is  usually  impossible  to  solve  exactly,  the  method  of 
moments  is  employed  to  generate  a  finite  dimensional  matrix  approximation  to 
Equation  (6.1).  The  process  begins  with  the  expansion  of  J  in  terms  of  basis 
functions,  in  this  case,  pulse  functions  that  are  defined  as 


P  (x,y) 
n 


(x,y)  e  cel  1  n 
(x,y)  t  cell  n 


(6.7) 


We  make  the  assumption  that  the  permittivity  of  each  cell  of  the  model  is 
constant,  and  number  the  cells  from  1  to  N.  Because  there  are  two  polarizations 
to  consider,  there  are  2N  unknowns  to  be  determined.  These  are  related  to  J  by 


J(x,y) 

"j^O 


L 

n=  1 


{x  X 


On(x,y) 


y  Yn  Pn(x,y)} 


(6.8) 


where  X  and  Y  are  unknown  coefficients.  From  Equation  (6.8),  we  obtain 
n  n 


curl  curl  A 


JUJE 


x  Y 


0 


92f 

n 

n  3x3y 


32  f 

<  —I 

n  ,  2 

3y 


32  f 

*  ’K  1737  -  Yn 


32  f 

_ £11 

2  J  1 
3x 


(6.9) 


where 


f  ( x ,y )  =  JJ  ~  H^2)(kR)  dx'dy'  (6.10) 

"  cell  0 

n 


Equation  (6.1)  may  be  separated  into  x  and  y  vector  components  and  enforced  at 
the  center  of  each  of  the  cells  to  yield  the  2Nx2N  system 
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,inc 


x  m  e 


32  f 


32  f 


x.  -  \  K  *  l  X"  TT 

n=l  3  n=l  3y 


m  =  1  ,2  ....  ,N 


N 

l  X 

“  n 


32f 


l  \ 


32f 


(6.11) 


Jy  m  £  -  1  m  L.  n  3x3y  ,  n  »  2 

3  rm  n=l  n=l  3x 


m  =  1 ,2,  .  .  .  ,N 


An  essential  feature  of  the  above  formulation  is  that  fn(x,y)  can  be  approxi¬ 
mated  analytically  and  then  differentiated  to  yield 


r. 


-> 

3*"  f 


3x 


■=—  (x  ,y  ) 
z  mm 


. 


ilti  H.  ( ka) 
4  1 


(x  ,y  )  £  cell  n 
m  -'ra 


juaJ^(ka) 


LH0(kp)  ~  +  Hj(kp)fZ — (xm,ym)  t  cell 
P 


2  2 
-  x 


(6.12) 


r 

2  0 
3*  f 

(xm,ym)  =/  j^aJ^ka)  f 
2 

V 


(x  y  )  £  cell  n 
tn  J  na 


[HQ(kp)  *  H}  ( kp)  )  ]  (xm.ym)  t  tell  n 


(6.13) 


32f 


32f 


— —  may  be  found  from  - by  exchanging  x  and  y  in  the  above  expression.  If 

3y  3x2  •  t 

the  suppressed  time  dependence  is  eJ  ,  Hq  and  represent  the  zero  and  first- 
order  Hankel  functions  of  the  second  kind.  An  assumption  involved  in  the  above 
process  is  that  the  cells  are  able  to  be  approximated  by  a  circular  cross  sec¬ 
tion  of  radius  'a.'  The  variables  (x  ,y  )  refer  to  the  center  of  cell  m, 

_  m  m 

n  i 

x=x  -x  ,  y-y  -y  ,  and  p  *  /x  +  y  . 
m  n  m  n 

6.3.  The  Symmetry  Structure  Imposed  by  a  Lattice  Geometry 

The  above  expressions  can  be  used  with  any  configuration  of  cells  in  a  model. 
For  geometries  of  the  type  depicted  in  Figure  6.1,  the  cells  are  arranged  in  a 
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two-dimensional  lattice.  Upon  investigation,  Equations  (6.12)  and  (6.13)  are 
seen  to  preserve  the  translational  structure  of  the  lattice,  so  that  if  the  num¬ 
bering  system  were  changed  to  the  type  depicted  in  Figure  6.2, 


j2fn  r  .  a2f22  r  . 

3x2  2'y2  '  3x2  3'y3 


(6.14) 


Thus,  the  efficient  numbering  system  of  Figure  6.2  enables  us  to  exploit  the 
two-dimensional  structure.  If  NX  and  NY  are  used  to  denote  the  lattice  dimen¬ 
sions,  Equation  (6.11)  can  be  expressed  in  the  form 


>  j6  3, 

.  .„-l; 


NY-1  NX-1 

I  I  v 

m=0  n=0 


GYX  .  + 

nm  a-n , B-m 


NY-1  NX-1 

l  l  x, 

m=0  n=0 


GXX 

nm  a-n , B-n 


'r  aB 
.  .0-1 


NY-1  NX-1 

I  V 

m=0  n^O 


X  GYX„  „  m  + 
nm  a-n , B-m 


NY-1  NX-1 

I  I  V. 

m=0  n=0 


GYY 

nm  a-n ,  8-ra 


a  =  0,1,..., NX- 1 


(6.15) 


6  =  0,1,.. . ,NY-1 

which  clearly  illustrates  the  symmetries  present  in  the  discrete  system. 
Equation  (6.15)  is  equivalent  to  the  2Nx2N  matrix  equation  (where  N=NXNY) 


(6.16) 


Each  of  the  blocks  A,  B,  and  D  is  itself  an  NX-th  order  block-Toepl itz  or 
almost  block-Toepl itz  matrix,  whose  individual  entries  are  Toeplitz  (or  almost 


Toeplitz)  matrices  of  order  NY.  (The  system  could  be  configured  so  that  the 
block- Toepl itz  systems  are  order  NY  and  their  entries  order  NX  Toeplitz 
systems.)  For  instance,  A  may  have  the  form 
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A  = 


"  (1) 

§0 


Si 


S2 


=NX 


Si 


,(2) 


Si 


S2 


So 


l 

(3) 


Snx 


(Nx) 


where  each  of  the  off-diagonal  entries  are  Toeplitz  matrices 


gg  gj  go  •  •  •  SfJY 

gl  g0  gl 

g2  gj  gg 


SNY 


80 


If  the  original  dielectric  cylinder  is  inhomogeneous,  i.e.,  if  s  varies 
cell  to  cell,  the  diagonal  entries  of  A  ate  not  exactly  Toeplitz,  and  are 
referred  to  as  almost-Toepl  itz  . 


(6.17) 


(6.18) 


from 
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As  compared  to  the  system  as  written  in  Equation  (6.11),  the  significant 
redundancy  features  of  the  latter  matrix  equation  permit  the  fully  populated 
2Nx2N  matrix  to  be  generated  from  only  51  entries.  Thus,  the  equation  is  an 
excellent  candidate  for  iterative  solution. 

6.4.  Fast-Fourier  Transform  Implementation  of  the  Matrix  Operator 

An  iterative  solution  of  the  above  matrix  equation  can  be  obtained  using  one 

of  the  algorithms  of  Chapter  2.  These  algorithms  do  not  require  the  matrix 

operator  explicitly,  only  that  it  (and  its  adjoint  or  approximate  inverse)  be 

available  to  operate  on  given  vectors  throughout  the  iterative  process.  Thus, 

the  manner  in  which  the  operator  is  implemented  is  of  no  consequence  to  the 

algorithm,  and  typically  is  accomplished  in  as  storage-efficient  manner  as  is 

possible.  One  approach  is  simply  to  perform  the  matrix  mul t ipl icat ion  directly, 

2 

which  requires  N  complex  operations  and  a  storage  of  2.5N  complex  variables, 
where  N=2(NX  *  NY). 

Since  the  matrix  operator  is  primarily  a  superposition  of  discrete  con¬ 
volution  operations,  it  can  be  implemented  with  the  aid  of  the  fast-Fourier 
transform  (FFT)  algorithm,  as  described  in  many  texts  { 4 2 1  ,  (43)  and  in  Chapter 
5.  For  the  two  dimensional  non-periodic  geometry  involved  in  this  example,  the 
task  requires  each  of  the  four  discrete  convolution  operations  to  be  replaced  bv 
three  successive  FFT  operations,  with  an  additional  vector  multiply  operation 
required  in  each  case.  (Actually,  only  two  FFT  operations  are  necessary  per 
convolution  during  each  iteration,  as  one  could  be  performed  initially  and  the 
result  stored.  The  discussion  below  assumes  that  this  is  the  case).  The  advan¬ 
tage  of  using  the  FFT  is  the  computational  efficiency  for  large  systems,  since 
approximately  Mlog(M)  operations  are  necessary  to  compute  a  one-dimensional  FFT 
or  inverse  FFT  of  length  M,  where  the  logarithm  is  of  base  2.  Thus,  each  of  the 
two-dimensional  discrete  convolution  operations  contained  in  the  matrix  operator 
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will  require  approximately  2 (2NX- 1 ) ( 2NY- 1 ) log( 2NX- 1 ) ( 2NY- 1 )  +  ( 2NX- 1 ) ( 2NY- 1 ) 
operations  via  the  FFT,  as  compared  to  (NXNYNXNY)  if  the  convolution  is  per¬ 
formed  by  direct  multiplication.  (These  estimates  assume  that  the  two  dimen¬ 
sional  FFT's  are  obtained  from  one  dimensional  FFT's,  as  is  typical  practice. 
Furthermore,  the  computational  effects  of  zero-padding  are  included  in  the 
estimates . ) 

As  illustrated  in  the  above  discussion,  the  effects  of  zero-padding  impact 
the  computational  efficiency  as  well  as  the  storage  efficiency,  and  specific 
calculations  are  considered  in  order  to  judge  the  overall  efficiency  of  using 
the  FFT  algorithm  to  implement  this  particular  operator.  For  instance,  if  the 
model  in  question  is  a  10x10  Lattice,  direct  matrix  multiplication  requires 
40,000  operations  and  a  storage  of  500  elements.  The  FFT  assisted  implementation 
of  the  matrix  operator  requires  25,980  operations,  and  a  storage  of  more  than 
2166  elements.  Thus,  for  large  problems,  the  FFT  approach  is  computationally 
efficient,  but  always  requires  mc'-e  storage  than  the  direct  approach.  In  addi¬ 
tion,  the  overhead  imposed  by  the  FFT  subroutine  adds  to  the  overall  program 
length  and  the  available  storage  space,  which  have  not  been  incorporated  into 
tne  above  estimates. 

In  summary,  the  matrix  operator  can  be  implemented  either  as  a  direct  matrix 
multiplication  or  as  an  FFT-assisted  computation.  Depending  on  the  specific 
geometry  involved,  one  approach  may  be  more  efficient  than  the  other.  In  a  given 
situation,  specific  parameters  describing  the  FFT  algorithms  available  for  use 
will  determine  the  actual  relative  efficiency.  It  may  be  worth  noting  that  the 
commercial  "array  processors"  frequently  supply  special  FFT  algorithms,  3nd 
their  availability  may  shift  the  relative  efficiency  toward  the  FFT-assisted 
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6.5.  Performance  of  the  Numerical  Approach 

The  performance  of  the  above  approach  can  be  studied  by  comparing  numerical 
solutions  and  exact  analytical  results  for  homogeneous,  circular  cylinders.  For 
instance,  three  models  of  circular  cross  sections  are  illustrated  in  Figure  6.3. 
Figures  6.4  to  6.12  show  values  of  the  total  electric  field  produced  by  the 
moment  method  approach  when  applied  to  these  models  and  compared  to  exact  solu¬ 
tions  based  on  an  eigenfunction  expansion  [94,.  In  this  case,  the  cylinder  under 
consideration  is  homogeneous  and  lossy,  with  £^=2 . 56- j2 . 56 ,  and  has  a  circum¬ 
ference  of  two  free-space  wavelengths.  As  the  models  are  refined  with  smaller 
cell  sizes,  the  numerical  solutions  appear  to  approach  the  exact. 

The  computer  program  that  generated  these  data  used  a  mixed  radix  FFT 
algorithm  to  implement  the  matrix  operator,  and  required  3.4,  24.2,  and  72.4  CPU 
seconds  of  execution  time  on  the  CDC  Cyber  175  computer  to  solve  the  N=42, 

N=202 ,  and  N=512  systems  using  the  CGM  algorithm  of  Chapter  2.  For  large  N,  this 
approach  clearly  exceeds  the  efficiency  of  the  MER  approach  discussed  in  Chapter 
-t ,  however,  execution  time  data  -an  be  misleading  if  not  interpreted  correctly. 
The  N=  312  equation  only  required  18  iterations  to  converge,  and  used  approxima¬ 
tely  3.4  seconds  of  CPU  time  per  iteration.  The  fast  convergence  is  related  to 
the  large  cell  density  in  use,  as  explained  in  Chapter  3. 

For  values  of  i  (  exceeding  10,  the  accuracy  of  the  numerical  solutions  is 
poor,  regardless  of  cell  density.  Furthermore,  the  convergence  of  the  CGM 
algorithm  becomes  less  and  less  rapid  as  j  e  j  increases,  suggesting  that  the 
system  becomes  ill-conditioned  for  large  |  ►:  |  .  This  seems  to  be  true  whether 
is  comp  1 ex-va 1 ued  (representing  a  lossy  media)  or  not.  For  £r=76- j2 78 ,  a  value 
often  used  to  model  biological  media  (81  j,  the  approach  fails  to  produce  numeri¬ 
cal  solutions  il.at  exhibit  any  kind  of  convergence  behavior  as  larger  cell  den- 
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Figure  b.b.  Comparison  of  the  exact  solution  and  the  numerical  solution 

obtained  with  the  256  cell  model.  The  cell  density  is  222  cells/A^. 
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Figure  6.7.  Comparison  of  the  exact  solution  and  the  numerical  solution 
obtained  with  the  21  cell  model. 
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Figure  6.9.  Comparison  of  Che  exact  solution  and  the  numerical  solution 
obtained  with  the  256  cell  model. 
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Fi^ur^  5.10.  Comparison  of  the  exact  solution  and  the  numerical 
obtained  with  the  21  cell  model. 
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sities  are  used.  The  reason  for  the  failure  of  the  approach  for  large  j  j  is 
not  as  yet  understood.  However,  it  is  not  a  direct  consequence  of  the  general 
DCMoM  procedure  or  the  iterative  solution.  It  may  be  a  result  of  the  integral 
equation  or  the  specific  basis  and  testing  functions  used  within  the 
d  iscret izat ion . 

6.6.  Summary 

The  DCMoM  procedure  can  be  applied  to  discretize  the  electric  field  integral 
equation  representing  scattering  by  inhomogeneous,  lossy  dielectric  cylinders. 
The  resulting  matrix  equation  is  convenient  to  solve  iteratively  because  of  the 
discrete-convolutional  redundancy  built  into  the  system.  In  this  manner,  effi¬ 
cient  use  can  be  made  of  available  fast-access  memory,  enabling  the  treatment  of 
electrically  larger  scatterers  than  possible  by  the  conventional  approach  [80] 
with  the  same  storage  constraints. 

From  a  comparison  of  numerical  solutions  and  exact  splutions  for  homogeneous, 
circular  cylinders,  the  process  appears  to  produce  accurate  solutions  for  |  | 

less  than  10.  The  numerical  solutions  for  problems  with  large  |e  |  fail  to  show 
convergence  behavior  as  higher  cell  densities  are  used  in  the  discretization 
procedure,  and  do  not  appear  to  be  reliable  approximations  to  the  desired  solu¬ 
tion.  The  reasons  for  the  failure  of  the  method  are  not  as  yet  understood. 
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7.  ITERATIVE  ANALYSIS  OF  HOLLOW,  FINITE-LENGTH  CIRCULAR 
CONDUCTING  CYLINDERS 

7.1.  Introduction 

Most  of  the  examples  presented  in  previous  chapters  involved  a  moment-method 
discretization  with  simple  puLse  basis  and  Dirac  delta  testing  functions.  In 
this  chapter,  a  DCMoM  discretization  using  higher-order  basis  and  testing  func¬ 
tions  is  applied  to  the  analysis  of  wave  scattering  by  a  finite- length ,  hollow, 
circular  conducting  cylinder. 

Because  circular  cylinders  are  periodic  in  the  azimuthal  variable,  all  quan¬ 
tities  can  be  represented  by  a  Fourier  series.  As  a  result,  Fourier  harmonics 
decouple  and  can  be  found  independently,  then  superimposed  to  produce  the  total 
solution.  This  procedure  is  known  as  a  "body  of  revolution"  formulation,  and  is 
desirable  because  the  matrix  equation  that  must  be  solved  for  each  of  the  har¬ 
monics  is  much  smaller  than  the  single  equation  arising  from  conventional  tech¬ 
niques  [95],  [96].  The  geometry  under  consideration  is  restricted  to  circular 
cylinders;  thus,  the  integral  equations  for  each  harmonic  can  be  discretized 
using  the  DCMoM.  Since  the  resulting  system  has  discrete-convolut  ional  sym¬ 
metries,  an  iterative  solution  could  require  significantly  less  storage  than  a 
conventional  direct  solution  and  add  to  the  savings  already  achieved  by  the  use 
of  the  "body  of  revolution"  formulation.  Therefore,  the  combination  of  the 
"body  of  revolution"  formulation  and  the  DCMoM  discretization  permits  the  treat¬ 
ment  of  significantly  larger  scatterers  than  otherwise  possible  with  the  same 
storage  constraints. 

A  different  approach  for  the  iterative  analysis  of  finite  circular  cylinders 
was  presented  by  Hurst  and  Mittra  (51 I.  Their  approach  was  based  upon  a 
straightforward  DCMoM  discretization  without  the  "body  of  revolution"  for¬ 
mulation  and  required  the  solution  of  one  large  matrix  equation  as  opposed  to  a 
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solution  of  a  number  of  smaller  order  systems.  An  example  involving  a  cylinder 
with  length  3.0  Aq  was  treated  via  the  iterative  solution  of  one  matrix  equation 
of  order  4032  [51].  Using  the  "body  of  revolution"  formulation,  the  same  problem 
could  be  treated  (to  the  same  degree  of  accuracy)  by  the  solution  of  a  number  of 
systems  of  order  63.  The  trade-off  in  computational  efficiency  between  these  two 
approaches  depends  on  the  number  of  Fourier  harmonics  excited  by  the  incident 
field  and  on  the  radius  of  the  cylinder.  In  general,  the  "body  of  revolution" 
approach  will  be  much  more  efficient  than  the  other  if  the  incident  waves  are 

propagating  in  a  direction  nearly  parallel  to  the  axis  of  the  cylinder,  for  few 

harmonics  are  excited  in  this  case. 

7.2.  Formulation  of  the  Matrix  Equation 

The  geometry  under  consideration  is  shown  in  Figure  7.1.  If  the  vector  com¬ 
ponents  of  the  current  density  are  expanded  in  Fourier  harmonics 
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If  the  cells  Ln  the  model  as  depicted  in  Figure  7.1  are  identical  in  size,  the 
discretization  process  reduces  to  a  form  of  the  DCMoM  and  produces  a  matrix 
equation  for  the  m-th  harmonic  of  the  form 
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The  sub-mac r ices  A  and  D  are  Toeplitz,  i.e. 
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where  the  superscript  "T"  denotes  matrix  transpose. 

The  computation  of  I™  requires  two-dimensional  numerical  integration,  and  is 
complicated  by  the  presence  of  a  singularity  in  the  integrand  when  q=0. 
Techniques  for  computing  1“  are  available  in  the  literature  |68),  [97]  .  For  the 
examples  to  follow,  numerical  integration  was  implemented  using  standard  adap¬ 
tive  library  routines. 

The  formulation  of  the  above  matrix  representation  is  identical  to  a  general 
approach  presented  by  Glisson  and  Wilton  ( 9  7 1  ,  if  their  approach  is  specialized 
to  the  geometry  of  Figure  7,1.  Detailed  information  about  the  matrix  elements, 
symmetries  in  the  equations,  and  the  decomposition  of  a  plane  wave  field  into 
Fourier  harmonics  may  be  found  in  the  literature  [96],  |97). 

7.3.  Incorporation  of  an  Impedance  Boundary  Condition 

The  above  integral  equation  is  based  upon  the  boundary  condition 


which  is  to  be  enforced  on  the  surface  of  a  perfect  electric  conductor  (PEC).  A 
generalization  of  this  condition  is 


E 

tan 


Z  J 
s 


(7.29) 


which  happens  to  coincide  with  the  classical  definition  of  surface  impedance  for 
an  imperfect  but  thick  conductor  [98].  A  boundary  condition  of  this  same  form 
has  also  been  used  under  certain  conditions  to  describe  electrically  thin 
materials,  and  may  apply  to  the  hollow  cylinder  geometry  under  consideration 
[99],  [100]. 

If  the  surface  impedance  Zg  is  invariant  to  azimuthal  variation,  boundary 
conditions  of  the  form  of  Equation  (7.29)  can  be  incorporated  into  the  above 
integral  equations  without  affecting  the  significant  symmetry  features  of  the 
DCMoM  matrix  equation.  The  effects  of  the  additional  terms  in  Equations  (7.3) 
and  (7.4)  are  to  modify  the  diagonal  entries  of  A  and  D,  and  may  destroy  the 
Toeplitz  natu.e  of  these  sub-matrices.  However,  the  off-diagonal  entries  retain 
the  discrete-convolutional  symmetries  which  permit  efficient  iterative  treatment 
of  the  system. 

7.4.  Performance  of  the  Method 

To  test  the  efficiency  of  the  above  approach,  two  computer  programs  were 
created  using  the  iterative  CGM  algorithm  to  solve  Equation  (7.13).  The  first 
program  used  explicit  matrix  multiplication  to  implement  the  operator  and 
adjoint  operator  required  by  the  CGM.  The  other  program  used  a  mixed-radix  FFT 
algorithm  to  perform  the  discrete  convolutions,  as  described  in  Chapters  5  and 
b.  As  expected,  the  FFT-assisted  approach  is  more  efficient  if  the  order  of  the 
matrix  equation  is  large.  However,  for  the  specific  computer  codes  in  use,  the 
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trade-off  between  the  two  approaches  was  not  as  favorable  as  expected.  For 
instance,  when  used  to  treat  a  99  x  99  system  of  equations  on  the  CDC  Cyber  175 
computer,  the  FFT-assisted  approach  required  approximately  0.35  second  of  com¬ 
putation  time  per  iteration,  as  opposed  to  0.24  second  per  iteration  for  direct 
matrix  multiplication.  For  a  399  x  399  system,  the  FFT-assisted  approach 
requires  about  1.8  seconds  per  iteration,  as  opposed  to  4.0  seconds  per  itera¬ 
tion  for  direct  matrix  multiplication.  (These  execution  times  should  be  com¬ 
pared  to  those  presented  in  Chapters  4  and  6.)  Although  the  efficiency  of  the 
FFT-assisted  approach  could  be  improved  by  using  a  radix-2  algorithm  [50],  the 
approach  is  obviously  not  as  efficient  as  expected  for  matrix  equations  of  rela¬ 
tively  small  order.  Furthermore,  use  of  the  FFT  algorithm  increases  storage 
requirements  somewhat.  Thus,  the  FFT-assisted  approach  may  not  be  beneficial 
unless  the  matrix  equations  of  interest  are  of  large  order  and  the  additional 
storage  constraints  do  not  exceed  the  limits  of  a  given  computer. 

Figure  7.3  shows  the  maximum  magnitudes  of  the  surface  current  density 
induced  on  a  large  perfectly  conducting  cylinder  by  an  axially  incident  plane 
wave,  according  to  the  above  integral  equations.  This  result  required  the  solu¬ 
tion  of  a  single  399  *  399  system,  and  the  CGM  algorithm  produced  a  satisfactory 
solution  after  79  iterations.  The  current  density  shows  interference  effects 
caused  by  superimposing  interior  and  exterior  currents,  as  is  necessary  when 
using  the  electric  field  integral  equation  to  model  thin  structures.  There  is 
no  known  exact  analytical  solution  for  the  finite  cylinder  geometry,  and  thus 
the  accuracy  of  the  numerical  result  is  not  determined. 

7.5.  Summary 

An  electric  field  integral  equation  representing  scattering  from  hollow, 
f  inite- length  circular  cylinders  can  be  discretized  using  the  DCMoM  procedure 
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combined  with  a  "body  of  revolution"  formulation.  The  system  for  each  Fourier 
harmonic  is  readily  solved  iteratively,  so  that  electrically  large  cylinders  can 
be  treated  with  minimal  storage  requirements. 
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8.  SUMMARY  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

This  report  describes  the  application  of  iterative  algorithms  to  the 
numerical  soLution  of  electromagnetic  scattering  problems.  Iterative  algorithms 
can  sometimes  be  incorporated  into  conventional  techniques  in  such  a  manner  as 
to  save  considerable  computer  storage  and  permit  the  treatment  of  larger 
discrete  systems  within  the  same  storage  constraints.  Issues  addressed  in  con¬ 
nection  with  iterative  methods  include  the  iterative  algorithms  available,  their 
performance  in  theory  and  in  practice,  and  the  various  ways  in  which  they  can  be 
incorporated  into  numerical  analysis.  Three  iterative  algorithms  related  to 
the  conjugate  gradient  method  are  presented  in  Chapter  2,  and  the  performance  of 
these  algorithms  when  applied  to  typical  electromagnetic  scattering  problems  is 
described  in  Chapter  3.  One  way  of  implementing  iterative  algorithms  is  the 
matrix-element  regeneration  (MER)  approach  discussed  in  Chapter  4.  The  MER 
requires  no  special  symmetries  in  the  matrix  equation,  but  is  not  as  efficient 
as  approaches  based  on  symmetries.  Chapter  5  discusses  two  discretization  pro¬ 
cedures  for  building  symmetries  into  the  matrix  equations,  and  detailed  examples 
of  this  type  of  approach  are  presented  in  Chapters  6  and  7. 

Since  the  goal  of  iterative  methods  is  numerical  efficiency,  execution  time 
data  have  been  emphasized  where  pertinent.  Whenever  possible,  numerical  solu¬ 
tions  have  ueen  compared  to  exact  analytical  results  for  verification,  and  their 
accuracy  was  studied  for  different  models  of  the  same  scatterer. 

Because  of  the  emphasis  on  discretizations  that  build  symmetry  into  the 
matrix  equations,  it  is  apparent  that  the  problems  best  suited  for  iterative 
solution  are  usually  those  involving  simple  geometric  shapes.  Many  additional 
cases  arise  in  which  parts  of  a  given  scatterer  will  conform  to  the  type  of 
geometry  easily  treated  iteratively,  and  other  parts  will  not.  For  these 
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problems,  iteration  may  be  efficient  if  a  large  enough  part  of  the  system  matrix 
can  be  made  to  possess  the  necessary  symmetry  features.  Sometimes,  it  may  be 
necessary  to  resort  to  rather  crude  modeling  in  order  to  achieve  this  type  of 
symmetry.  Depending  on  the  problem  under  consideration,  the  MER  approach  may  be 
a  viable  alternative  to  crude  modeling  techniques. 

Although  the  use  of  iterative  algorithms  based  on  gradient  methods  has 
eliminated  the  problems  with  diverging  numerical  solutions  that  plagued  earlier 
researchers,  faster  convergence  is  always  a  desired  goal.  In  practice,  the 
algorithms  described  in  Chapter  2  may  be  expensive  if  applied  to  large  systems. 
It  is  expected  that  future  efforts  in  the  fields  of  engineering  and  computer 
science  will  improve  the  available  algorithms,  and  improvements  could  and  should 
be  incorporated  into  existing  methods  as  quickly  as  they  become  available. 
Existing  algorithms  based  on  preconditioning  the  matrix  equation  require 
knowledge  concerning  the  eigenvalues  of  the  typical  systems  that  arise.  Based 
upon  the  performance  of  the  algorithms  as  illustrated  in  Chapter  3,  there  appear 
to  be  similarities  in  the  eigenvalue  behavior  of  different  scattering  problems 
that  could  be  identified  and  used  in  preconditioning  algorithms.  This  may  moti¬ 
vate  an  investigation  into  the  eigenvalue  'structure  of  the  typical  integral 
equations  of  electromagnetics. 

The  topic  of  multiple  right-hand  sides  remains  3  vexing  problem  for 
electromagnetic  scattering  problems,  and  deserves  additional  attention. 

However,  as  discussed  in  Chapter  3,  the  goals  of  fast  convergence  and  ability  to 
treat  multiple  right-hand  sides  are  at  odds,  and  algorithms  well-suited  to  the 
first  will  probably  be  poor  at  the  second.  It  appears  that  attempts  to  treat 
multiple  right-hand  sides  should  be  based  upon  orthogonal  expansions  with  more 
flexibility  than  those  produced  by  the  CGM,  which  are  geared  to  one  solution  of 


a  given  system. 
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The  state-of-the-art  techniques  for  the  solution  of  discrete  systems  need 
to  reflect  current  progress  in  hardware  as  well  as  in  software.  The  avail¬ 
ability  of  large  blocks  of  random  access  memory  masquerading  as  dummy  disks  and  the 
enhancement  of  computer  architecture  to  allow  larger  chunks  of  directly 
addressable  memory  both  have  an  immediate  impact  on  the  size  problem  which  can 
be  treated,  and  the  efficiency  of  direct  and  iterative  algorithms.  These  issues 
must  be  kept  in  mind  when  considering  a  method  for  a  specific  problem.  Of 
course,  it  is  expected  that  the  size  problem  that  one  would  like  to  be  able  to 
solve  will  always  exceed  the  capacity  of  existing  machines. 

It  is  evident  from  this  work  that  there  are  many  alternatives  and  trade¬ 
offs  to  be  addressed  when  considering  an  iterative  computational  method  for  the 
solution  of  electromagnetic  scattering  problems.  The  character  of  a  specific 
problem,  the  computational  facilities  available,  -and  the  desired  accuracy  in 
modeling  and  in  the  numerical  solution  will  determine  the  best  approach  for  a 
given  problem. 
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